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\Q 1 ABSTRACT 

O . 

We define three requirements for accurate simulations that attempt to model cir- 
cumstellar disks and the formation of collapsed objects (e.g. planets) within them. 
O !■ First, we define a resolution requirement based on the wavelength for neutral stabil- 

D ■ ity of self gravitating waves in the disk, where a Jeans analysis does not apply. For 

\ particle based or grid based simulations, this criterion takes the form, respectively, of 

OO . a minimum number of particles per critical ('Toomre') mass or maximum value of a 

'Toomre number', T — 6x/Xt, where the wavelength, At, is the wavelength for neu- 
tral stability for waves in disks. The requ irements are ana l ogues of t he conditions for 
cloud collapse simulations as discussed in iBate &: Burkerfl l)1997j) and iTruelove el al. I 
^ ' l)1997|) . where the required minimum resolution was shown to be twice the number of 

neighbors per Jeans mass or 4-5 times the local Jeans wavelength, Aj, for particle or 
grid simulations, respectively. 

We apply our criterion to particle simulations of disk evolution and find that in 
order to prevent numerically induced fragmentation of the disk, the Toomre mass must 
be resolved by a minimum of six times the average number of neighbor particles used. 
\ We investigate the origin of the apparent discrepancy between the number of particles 

^ — required by the cloud and disk fragmentation criteria and find that it is due largely 

■ to ambiguities in the definition of the Jeans mass, as used by different authors. We 

reconcile the various definitions, and when an identical definition of the Jeans mass is 
Q ■ used, the condition that </< 1/4 in the Truelove condition is equivalent to requiring 

\ about 10-12 times the average number of neighbor particles per Jeans mass in an SPH 

simulation, reducing the difference between simulations of disks and clouds to about 
two. While the numbers of particles per critical mass are similar for both the Jeans 
and Toomre formalisms, the Toomre requirement is more restrictive than the Jeans 
requirement when the local value of the Toomre stability parameter Q falls below 
about one half. 

^ . Second, we require that particle based simulations with self gravity use a variable 

gravitational softening, in order to avoid inducing fragmentation by an inappropriate 
choice of softening length. We show that using a fixed gravitational softening length 
for all particles can lead either to artificial suppression or enhancement of structure 
(including fragmentation) in a given disk, or both in different locations of the same 
disk, depending on the value chosen for the softening length. Unphysical behavior can 
occur whether or not the system is properly resolved by the new Toomre criterion. 

Third, we require that three dimensional SPH simulations resolve the vertical 
structure with at least ~ 4 particle smoothing lengths per scale height at the disk 
midplane, a value which implies a substantially larger number per vertical column 
because the disk itself extends over many scale heights. We suggest that a similar cri- 
terion applies to grid based simulations. We demonstrate that failure to meet this cri- 
terion leads to underestimates in the midplane density of up to 30-50% at resolutions 
common in the literature. As a direct consequence, gas pressures will be dramatically 
underestimated and simulations of self gravitating systems may artificially and erro- 
neously inflate the likelihood of fragmentation. We outline an additional condition on 
the vertical resolution in simulations that include radiative transfer in order to ensure 
a correct description of the cooling, specifically that the temperature structure near 
the disk photosphere must be well resolved. As an example, we demonstrate that for 
an isentropic vertical structure, the criterion translates to resolution comparable to 

, — . _„__ „ . „ H/20 near the disk photosphere, to avoid serious errors in transfer rates of thermal 
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energy m and out of the disk. 

Finally, we discuss results in the literature that purport to form collapsed objects 

and conclude that many are likely to have violated one or more of our criteria, and 

have therefore made incorrect conclusions regarding the likelihood for fragmentation 
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1 INTRODUCTION 

The formation of collapsed objects, both in the context of 
the collapse of molecular cloud cores and in the later context 
of the clumping of material in a circumstellar disk is very 
difficult to model numerically because of the huge range of 
spatial scales involved. Accurate simulation of the collapse 
usually requires that all of these scales be well resolved, if 
the result is not to be contaminated by numerically induced 
fragmentation. 

F or three dimensional (3D) simulations. iTruelove et al. I 
il997T) defined a minimum resolution condition for the nu- 
merical validity of a simulation that models the collapse of a 
molecular cloud core usin g a grid based hy d rodyn amic code. 
Contemporary work by iBate fc Burkertl <ll997l . hereafter 
BB97) has examined necessary resolution conditions for the 
collapse of a similar system in the context of Smoothed Par- 
ticle Hydrodynamics (SPH) simulations, and also discussed 
the influences that choices of the form of gravitational soft- 
ening may have on the results. Both of these works define 
minimum resolution criteria in the context of a Jeans col- 
lapse o f a gas cloud resolved in 3D, but while ITruelove et al. I 
il997fl observe that fragmentation is enhanced by a failure 
of the criterion, BB97 only observe enhanced fragmentation 
if the gravitational softening used in their particle based 
simulations is smaller than the hydrodynamic smoothing. 
Otherwise, they observe that fragmentation can be unphys- 
ically d elayed in a system where it is known to occur. Later 
work of Hubbcr . Goodwin fc W hitworth ( 2006) explores the 
problem of Jeans collapse in isolation. They find that SPH 
simulations of systems with initial conditions like that of the 
original linearized analysis do not fragment artificially, even 
when under resolved. 



A number of recent wor ks in the fie l d of planet 



and brown dwarf formation JNelson et al. 1 Il998l. 


2000 


Pickett et al. 1 ll99Sl. l2000al. 120031 iBoss! Il998l. 


200C 


2001 120041: iMaver et al. 1 12001 12004 


Rice et al. 1 


2003 


Armitaee & Hansen! Il999t iLufkin et al. 


2004) have dis- 



cussed the formation of planets and brown dwarfs via 
gravitational fragmentation in disks. In disks systems like 
those discussed, the Jeans formalism used to develop the 
ITruelove et al. I Il997l and BB97 criteria is not valid both 
because the disk scale height is usually small compared to 
the Jeans wavelength and because of the existence of shear, 
which plays a role as important as self gravity for the grow- 
ing structures. As yet however, no analogous resolution cri- 
terion exists for disks, which may be modeled in either fully 
in two dimensions (2D), effectively integrating over the ver- 
tical coordinate, or in 3D, for which no vertical integration is 
assumed, but for which the Jeans wavelength based criterion 
still may not apply. 

In addition to requirements that simulations resolve the 
wavelengths of the relevant instabilities sufficiently, simula- 
tions must also satisfy a number of other criteria if they 
are to be believed. Approximations made outside the realm 
of the physical model, perhaps used to model the behav- 
ior of unresolved or poorly resolved phenomena, must not 
drive the results of the simulations themselves. In the con- 
text of simulations of disks, particularly those using parti- 
cle based hydrodynamic methods such as SPH, an impor- 
| Current Address: Los Alamos National Laboratory, X-2 MS tant consideration will be the implementation of gravita- 

T087, Los Alamos NM, 87545, USA tional softening. At a more fundamental level, the algorithm 
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used to solve the equations used to model physical system 
must do so accurately, without becoming unstable or gen- 
crating large errors through some other sort of deficiency. 
A quick glance through the literature (see e . g. textbooks o f 
iHocknev fc Eastwoodfll988t iFletched Il997t iLeveaud I2OO2T) 
will demonstrate to even a casual observer that the study 
of numerical methods in relation to their stability is one 
in which extensive studies on many topics have been per- 
formed. Verification of methods on physically relevant test 
problems however, is comparatively more widespread but 
may often suffer from insufficient detail or generality. While 
many studies discuss the fidelity of the numerical solution on 
some simple or contrived test problem, often no test prob- 
lems sufficiently similar to the physical system under study 
can even be devised. 

In sect ion we first extend the previous 3D work of 
True love et al. and BB97 to self gravitating thin disk 

systems, then outline alternatives for gravitational softening 
in particle simulations, and the possible consequences each 
choice may have on results. In section |H] we define a test 
problem with which we determine appropriate values for 
resolution for particle simulations. We continue in section 
2] with a discussion of the gravitational softening assumed 
in particle simulations, and the specific numerical issues en- 
countered in a study of disks. Next, in section|SJ we discuss 
the application of the criterion to simulations done in 3D, 
and define a test problem to validate the accuracy of numer- 
ical codes attempting to model disks in 3D. Using this test 
problem, we demonstrate that failure to resolve the vertical 
structure of disks will lead to large errors in the numeri- 
cal solution for the evolution of the entire physical system. 
Finally, in section |SJ we summarize our results and discuss 
them in the context of the models presented in the literature. 



2 NUMERICAL FACTORS AFFECTING THE 
RESULTS OF DISK SIMULATIONS 

A numerical simulation of any physical system may suffer 
from inaccuracies from several different origins. For example, 
an incorrect physical model or incorrect initial conditions 
will produce results irrelevant to the system being modeled. 
On the side of numerics, a shortcoming in the numerical 
method may erroneously trigger some physical process to 
become active in the evolution, where in reality, no such 
physical process is important. A shortcoming in the numeri- 
cal method may also trigger effects of purely artificial origin. 
In this section, we first describe the conditions under which 
we may expect a numerically induced, but physically based 
instability leading to fragmentation to be present in simu- 
lations of disks. Secondly, we discuss alternative treatments 
of the hydrodynamics and gravity in particle simulations on 
the spatial scales of the smoothing and softening lengths in 
particle simulations. We describe conditions under which we 
may expect their implementations to influence a simulation, 
possibly also leading to artificially induced fragmentation. 
Finally, we discuss the strengths and weaknesses of 2D and 
3D treatments of a system, the meaning of physical quanti- 
ties realized in a 2D approximation and consequences that 
may arise when the assumptions underlying one or the other 
treatment break down. 



2.1 Resolution criteria in cloud and disk 
fragmentation simulations 

A condition on the minimum resolution to ensure the col- 
lapse of a clou d is of physical rather than numerical origin 
was defined bv lTruelove et al. I Jl997ft . using the ratio of the 
local grid resolution, Sx, and local Jeans wavelength, Aj, in 
the fluid: 



J : 



6x 



(1) 



where 8x is the size of a grid cell and Aj is the local Jeans 
wavelength: 

/ 2\V2 

• (2) 

and c s is the sound speed, p is the volume density and G 
is the gravitational constant. To obtain a form more useful 
in particle based numerical methods (e.g. SPH) where the 
resolution element is a unit of mass rather than of length, 
BB97 us ed the Jeans mass as defined from energy consid- 
erations <lTohlinelll982t) for a homogeneous sphere, to define 
an analogous criterion for the maximum resolvable density. 
Generalizing their result to a gas with d internal degrees of 
freedom we find: 



Al 



Energy 
J 



_3_ 

4tt 



-, 1/2 



(G 3 p) 1/2 



(3) 



where the superscript 'Energy' is included to distinguish this 
definition of Mj from two others defined in the next section, 
7 is the ratio of specific heats and d is the number of de- 
grees of freedom, equal to 3 for a monotonic ideal gas. This 
equation yields a maximum resolvable density : 



Pmax 



_3_ / 5d V 
4tt \6jG ) 



(m p N reso ) 2 



(4) 



SPH 



where we equate the Jeans mass Mj to a sum of N TI 
particle masses, m p , that are required to resolve it. 

An exactly analogous stability condition can be made 
for rotationally supported (i.e. disk) systems using the local 
Toomre wavelength, At, 

T=£, (5) 
at 

where At is the wavelength which defines neutral stability 
in disks. We can derive At from the d ispersion re lation for 
waves in disks, whose solution (see e.g.. lLin fc Laull979l) has 
four branches: 



±k (l ± v/f-Q^I-^ 2 )) 



(6) 



corresponding to leading and trailing, short and long wave- 
length spiral density waves. The variables ko, Q and v are 
defined by 



fco 



vrGE 



ttGE 



and 



v = (^~ mQ ) (7) 



respectively. Physically, these variables are the wave num- 
ber, the well known Toomre Q parameter and the Doppler 

1 Note that t his equation in our previous conference proceeding 
iNelsonl J2003fl was erroneously stated and should be disregarded. 
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shifted pattern frequency of wave of symmetry m (i.e. with 
m spiral arms), normalized to the local epicyclic frequency, 
k, in the disk. They depend on the disk's surface density 
E, its orbit frequency, u, as well as the pattern frequency, Q 
and symmetry, m, (i.e. the number of spiral arms) of the spi- 
ral waves. Neutral stability is defined by the condition that 
the term under the square root is zero, for which the wave 
number is jfcj = fco. The critical wavelength corresponding 
to this wave number is: 



At — 



2cj_ 
GE ' 



(8) 



iBinnev fc Tremaind il987f) derive an expression for the 
longest wavelength that will be unstable in a given disk as 



Ac 



4vr 2 GE 



(9) 



however this wavelength is of limited relevance to the ques- 
tion of numerical resolution because it is neither the most 
unstable nor the shortest unstable wavelength in the disk. 
As noted by them, when the disk first becomes unstable 
(i.e. Q — 1), A c is a factor of two longer than At- Due to 
our interest in determining a critical resolution requirement, 
which requires the shorter length scale be used, we use the 
wavelength definition of eq. |5] rather than that of eq. |§] in 
defining our criterion. 

As in the 3D case, a form useful for particle based sim- 
ulations can be obtained, this time by defining a Toomre 
mass. Unlike the 3D case however, due to the more compli- 
cated geometry of the mass distribution, no easily derivable 
form of the Toomre mass may be determined from energy 
considerations analogous to that in eq. [3] We therefore use 
a circular volume element to define 



M T = 



G 2 E ' 



(10) 



A maximum resolvable surface density follows directly as: 



^max — 



VGV m p N It 



(11) 



where the symbols have the same meanings as in equation 2] 
above. Determination of appropriate values for T and A^eso 
for disk simulations will be the subject of section [3] 

It will frequently be useful to determine the resolution 
required for the same overall morphology but with vary- 
ing stability, as required for a parameter study varying a 
disk's Toomre Q value. Given the same overall morphol- 
ogy, Q becomes a function of sound speed only, which is 
also the only term in equation 1111 that will vary. A useful 
re-parameterization of this expression to illustrate the sen- 
sitivity of the resolution to the disk stability directly will 
therefore be to replace the sound speed with Q through its 
definition: 



tt 5 G 2 E 4 



Q 4 



m„Nr, 



(12) 



The maximum resolvable surface density is thus propor- 
tional to the fourth power of Q. In order to obtain iden- 
tical effective resolution, as quantified by identical values of 
E max , a simulation with Q = 1 must therefore be resolved 
with Q 4 « 5 more particles than a simulation with Q = 1.5. 



2.2 Consistent application of the criteria across 
simulations of different types 

The criteria in equations and 0] can easily be shown to be 
equivalent up to constant factors by casting the former as a 
Jeans mass. The mass inside a sphere of radius Aj/2 is: 



A I 



Sphere 
J 



47T /A 



24 (G 3 p)V2- ( 13 ) 

so that comparison of equations|3](determined from a similar 
spherical volume assumption) and equation 1131 yields the 
desired identifi cation. Yet another def inition of the Jeans 
mass is given m iKrumholz el al. I i2004Tl using a cubic rather 
than spherical volume element as: 



M Cube = p> 3 



3/2 



(G 3 p) 1/2 



(14) 



In each of these definitions, the constant factors differ. 
Specifically, the three forms Mf nergv : M s , phere : Mf ube , 
yield values of the Jeans mass in a ratio of 0.89 : 2.92 : 
5.57 relative to each other, assuming a monatomic ideal gas. 
Analogously for the 2D case, using a square area element, Ay 
rather than circular, would result in a definition of Toomre 
mass or density that is a factor 4/-7T larger than stated in 
equations 1101 and II II respectively. 

Although these constant factors may indeed be regarded 
as insignificant in a global sense since the basic functional 
dependence is the same, their quantification is important 
because it allows us to interpret the results of simulations 
obt ained using various ver sions of the criterion. For exam- 
ple, [^nielov^e^II] lll997l) found empirically that J<l/4 
was sufficient to suppress numerical instabilities in a test 
problem, which implies a mass per grid zone of at most 
m zone ~ J 3 Mj ~ Mj/64. For particle simulations, BB97 
showed that numerical instability could be suppressed in a 
similar test problem by resolving the local Jeans mass with 
particles less massive than m p ~ Mj/N Teso ss Mj/100. Al- 
though ostensibly very similar, these criteria differ by a hid- 
den factor of ~ 5.57/0.89 ~ 6.2 from each other due to the 
differing proportionality factors, with the particle based cri- 
terion being more conservative (requiring more particles to 
obtain equivalent resolution for a given mass distribution). 

In other words, the ~ 100 particle (w x2iV ne igh) con- 
dition of BB97 translates to a ~ 300 particle (~ x6iV no i g h) 
condition condition if we use the Jeans mass defined by equa- 
tion [HI or a ~ 600 particle (~ xl2A ne i g h) condition if we 
use the Jeans mass definition of equation 1141 In general, 
for equivalent resolution of the physical parameters of the 
flow, a particle simulation will require roughly ten times as 
many fluid elements a s a grid code like the one discussed in 
iTruelove et al. I dl997fl . With the same number of fluid ele- 
ments (particles or grid cells) a grid simulation will be able 
to resolve higher mass densities before becoming numerically 
unstable than a particle simulation. 



2.3 The relative sizes of the Jeans and Toomre 
wavelengths 

The critical wavelength for disks from equation |H| is linearly 
dependent on the temperature to mass ratio through the 
sound speed and mass density, while the Jeans wavelength 
in equation|5|is dependent only on its square root. Thus, the 
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Toomre criterion will be more strict on smaller spatial scales 
than the Jeans criterion, but less strict on larger scales. 

We can determine the relative sizes of the two wave- 
lengths, and the crossover point at which both are equal, 
from their ratio: 



At 

17 



PCs 



ttGE 5 



1/2 



(15) 



If we assume that the disk is near Keplerian so that w «, 
that the volume density and the surface density are related 
by p — E/2/_ff, where / ~ 1 is a coefficient specifying the 
exact proportionality between surface and volume densities, 
and that the local disk scale height is H — c s /Q (but see 
section IB~21 below') . then we can combine equation [5] and [H] 
to produce 



A 3 



(16) 



Thus if the local value of Q falls below about one half, as it 
may in regions beginning to fragment, the Toomre instability 
wavelength will be a more restrictive criteri on for numerical 
simula tions than the Jeans criterion used bv lTruelove et al~l 
(1997). 

On first inspection, eauation ll6l and the wavelength pro- 
portionalities that went into it would seem to be backwards. 
When the collapse is well underway (Q « 1), rotation 
ceases to matter so that the collapse should proceed ac- 
cording to a Jeans prescription. In other words, we would 
expect that the Jeans wavelength should be smaller than 
the Toomre wavelength. Indeed, this would be the case for 
a collapsing region that is much smaller than the disk scale 
height, but this equation shows that for the initial stages of 
collapse relevant to the analytic wavelength derivations, it 
is not the presence or absence of rotation that is relevant, 
but rather the dimensionality of the problem that plays the 
most important role. 



2.4 The application and applicability of the Jeans 
and Toomre criteria to numerical simulations 

The mathematical analysis leading to the Jeans resolution 
criterion is fully three dimensional, while that leading to the 
Toomre criterion is limited to two dimensions: the equations 
are integrated over the z coordinate. The difference is im- 
portant because the Jeans analysis may not be valid in a 
disk simulation (even those evolved using fully 3D models) 
because it assumes a homogeneous medium which is infinite 
in all three spatial coordinates. By definition, a disk struc- 
ture will violate this assumption since the matter will always 
be condensed into a midplane above and below which rela- 
tively little matter lies. The violation will be important if the 
Jeans wavelength is comparable to or larger than the disk 
scale height, determined for the local conditions in a given 
disk. In this case, even the initial conditions of a system may 
not satisfy the underlying assumptions of the analysis. Tak- 
ing the specific example shown in figure l2Ub of section lb. 21 
below, we find that Jeans wavelength is long compared to 
the disk scale height everywhere, so that the Jeans analysis 
is indeed inapplicable. 

On the other hand, the analysis leading to the Toomre 
wavelength is valid in the limit of a thin, rotating system for 



which a 'surface density' is a meaningful concept, whether 
or not a given simulation is actually performed in two or 
three dimensions. By construction of the analysis itself, disks 
fall into this category, so at least we may construct initial 
conditions that satisfy the underlying model. As in the case 
of the Jeans wavelength, by examining figure !2Ub . we see 
that the Toomre wavelength is also long in comparison to 
the disk scale height. In this case, the large ratio means only 
that the analytic assumptions become more valid, not less. 
Violation of the model assumptions may still be important 
with wavelengths comparable to a disk scale height scale 
because of the neglect of vertical motion and structure in 
the analysis. 

Violation of the assumptions in the analytic derivations 
may also occur at times after simulations have evolved for 
some time because the analyses assume that variations of 
all quantities from their initial values are small, even if the 
initial conditions satisfy all requirements of the linearized 
analyses. For example, a density perturbation may be only 
slightly enhanced from the local background in a fragment- 
ing molecular cloud, while the background potential is char- 
acterized by a comparatively steep overall gradient due to 
some large structure nearby. A second example may be that 
fluid velocities are large while other perturbations are small 
because the initial state was seeded with some spectrum of 
turbulent velocity perturbations. Finally, numerical simula- 
tions may not solve the equations of hydrodynamics accu- 
rately, or may do so with only poor fidelity for some prob- 
lems. 

Unfortunately, we will find that in the most relevant 
range of parameter space for disk simulations, both the 
Jeans and Toomre wavelengths reach values comparable to 
the disk scale height and local perturbations reach ampli- 
tudes beyond those for which linearized analyses are strictly 
valid. As with the original 3D analysis of collapsing clouds, 
we therefore resort to numerical simulations to determine 
the required 'safety factor' (i.e. the values of N IOSO , T' and 
'J' above for disk or cloud collapse analyses) for which we 
can be reasonably confident that the basic features of that 
analysis are not violated, even though it may be applied in 
a region where its assumptions may be called into question. 

Although there is only one criterion for disk evolution 
and one for cloud collapse, there are effectively five imple- 
mentations of them as applied to disk systems that could 
be used under different circumstances. First the equations^ 
or might be used directly: a 2D simulation could use the 
directly available surface density, or a 3D simulation could 
use the directly available volume density. The criteria could 
also be used indirectly using the disk scale height to con- 
vert between volume density and surface density. Finally, 
a surface density could be obtained from a 3D simulation 
by directly integrating over the z coordinate, for use in the 
Toomre criterion. We will show in section 15.21 that using 
the approximate indirect forms (i.e. making the volume or 
surface density conversion using the disk the scale height) 
can yield seriously discrepant values of the stability wave- 
length, perhaps leading to erroneous conclusions regarding 
the veracity of a given simulation. 
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2.5 Hydrodynamic smoothing, gravitational 

softening, and our implementations of them 

In all numerical simulations, the modeler would like to re- 
solve the largest range of spatial scales possible, so that both 
the smooth and highly inhomogeneous regions are accurately 
modeled. In this regard, a limit on the resolution will be the 
scale on which the gravitational and hydrodynamic forces 
can be resolved. For grid based simulations, resolution of 
both the hydrodynamic features in the flow and gravita- 
tional forces will be related to the local dimensions of the 
grid (see e.g. iFrvxell. Miiller fc Arnettlll99lt iPickett et al\ 

For particle simulations, the limits will be related to two 
length scales, one for gravity and one for hydrodynamics, 
each defining the spatial extent of the particles in different 
ways. In order to make clear that they are distinct from each 
other, we will make a distinction between the terms used 
to refer to each. Specifically we will refer to gravitational 
'softening' lengths and hydrodynamic 'smoothing' lengths 
for particles, to describe each effect. 



2. 5. 1 Smoothing 

Fluid quantities in a particle based simulation using SPH 
are reconstructed from the positions of the particles at a 
given time. Contributions of particles are weighted accord- 
ing to a smoothing function (the 'kernel') and all contribu- 
tions are summed to define each of the hydrodyn amic quan- 
tities at the current location of each particle jMonaghanl 
1992). Quantities such as density are calculated using the 
kernel directly, while forces due to pressure are calculated 
using its derivative. I n each case our implementation fol- 
lows the discussion in iBenzl jl990t) and, in the work pre- 
sented here, we use the now standard spline kernel of 
iMonaehan fe Lattanzid il985fl given by: 



1 



W(nM = — i \{2-vf 
% I o 



if < v < 1 
if 1 < v < 2 
otherwise. ; 



(17) 



Here, v is the number of dimensions, v = nj/hij , ry = n — 
r j I ! hij = ( hi + hj ) /2 and a is the normalization with values 
of 10/(77r) and 1/tt in two or three dimensions respectively. 

We will also study the effect of an important mod ifica- 
tion o f that derivative described by iThomas fc Couchmanl 
( 1992), which takes the form: 



if < v < 2/3 



h» +1 

*3 



-1 

-3v + !« 2 if 2/3 < v < 1 



"!(2-» 




% ; f 7<;< 2 7' (18) 



f/U ( / , ,. 

otherwise ; 

and acts to reduce unphysic al particle clumping iSteinmetzl 
ll99rJ: iThacker et aTlhoOOh . by providing a net repulsive 
pressure force even at zero separation. The modification af- 
fects the region v < 2/3 where in its unmodified form, the 
derivative decreases monotonically to zero as v itself ap- 
proaches zero. Price (2005, personal communication) notes 
that this modification may cause changes in the effective 
sound speed determined from linear analyses, because nor- 
malization conditions on the kernel's first and second deriva- 
tives are not satisfied. The normalization conditions will be 
exact only in the limits where the smoothing lengths ap- 



proach zero, neighbor counts approach infinity and the dis- 
tribution of neighbors over the kernel volume allows an ac- 
curate correspondence between a summed and an integral 
form of the equations. Since none of these three conditions 
are realized in practice, it is difficult to evaluate what con- 
sequences may result in actual simulations. In section 13.21 
we investigate the differences between using the standard 
and modified kernel derivatives on the outcomes of our sim- 
ulations. Detailed analyses of the numerical stability and 
errors resulting from this modified formulation are beyond 
the scope of this work and we defer such discussions to the 
future. 



2.5.2 Softening 

An important subtlety of particle simulations is to employ 
a gravitational softening length that removes the singular- 
ity in the force obtained when two particles in a simula- 
tion approach each other too closely, effectively modifying 
Newton's law of gravitation. Such a modification is equiva- 
lent to defining a mass density distribution for the particles, 
so that they include an assumption of some spatial extent, 
rather than that they are point like objects. Modifying the 
gravitational force is acceptable in hydrodynamic and col- 
lisionless A-body simulations because particles do not in 
fact represent single particles in the underlying system, but 
only a statistical representation of the local distribution of 
gas or particles. The question we face is how to implement 
an appropriate gravitational softening length in our simula- 
tions which, on one hand, prevents unphysically large forces 
from developing between particles, and on the other, is small 
enough to allow small scale features to develop in the flow. 

Two alternatives for implementing softening that are 
common throughout the literature are first to assume each 
particle has the mass distribution of a Plummer sphere, so 
that the force law is modified to the form of a Plum mer force 
law llRomecll994l . ll997llAthanassoula et al. l 200Cf) . or to use 
a spline based kernel as discussed bv lBenzl il990l) . who inter- 
prets the smoothing kernel used to realize the hydrodynamic 
quantities as a mass distribution. Two variants are common 
in the latter case, first to use a softening that varies with 
the local conditions (usually chosen to be identical to the 
SPH smoothing length of each particle) or second, to fix the 
softening to a single constant value at the beginning of the 
simulation. Exa mples o f varia b le kernel softening ar e foun d 
in codes used by Ben j (Il99(t) : ISteinmetz fc Mullerl il993ft . 
BB97 and our own previous work, while the TreeSPH , GAD - 
GET and GASOLINE cod es of iHernauist fc Kap i 19891) : 
ISpringel et oTI feOOll) and IWadslev et al. I (I2003T) . respec- 
tively, use fixed softening. 

We implement the variable kernel softening variant us- 
ing same kernel used for the SPH quantities given by equa- 
tion !17l For gravitational softening, the masses used to com- 
pute the acceleration is modified from its Newtonian form 
according to the prescription that a source particle's mass is 
reduced from its true value by a factor proportional to the 
volume enclosed by a sphere whose radius is the separation 
between the particles: 



(19) 



2tt J^' j Wudu if n d = 2; 
4tt J ry Wu 2 du if n d = 3 

where m p is the mass of the particle from which a force 
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contribution is to be calculated and nj is the number of 
spatial dimensions. The sink particle, on which forces are 
calculated, is assumed to be a point mass, so that for two 
particles of mass, mi and mi, the force exerted by particle 
1 on particle 2 is 



—Grh\m2i\ 



(20) 



where ri2 = ri — is the separation between the parti- 
cles one and two. Except for a change of sign, this definition 
is manifestly invariant to exchanging the identities of the 
two particles, and therefore conserves momentum exactly. In 
3D, and as the separation between two particles decreases, 
the gravitational force between them will also decrease, ul- 
timately to zero at jrj = 0, because m depends on r 3 . 

In 2D, m is proportional only to r 2 and the force instead 
approaches a non-zero constant value as the separation de- 
creases to zero, rather than zero. The physical reason for the 
discrepancy is that while equation ED implies a 2D structure 
for the mass distribution, eauation l2Ul retains a 3D structure 
for the forces they cause. On the scale of the interparticle 
separation, the scale height of the disk not be negligible, so 
that in fact it is truly three dimensional. We are aware of 
only one treatment that attempts to acc ount for the con- 
flict between the 2D and 3D behaviors jKollerlliooi) . by 
assuming a vertical structure and numerically integrating 
the contributions to the forces on a point mass over the z 
coordinate. Thereafter Koller applies the derived correction 
factor to the forces, tabulated as a function of distance. We 
expect Roller's treatment will not be generally suitable for 
particle simulations because it requires different tables for 
different vertical structures and because, even with the mod- 
ification, he finds that a Plummer-like softening term is still 
required, though it can be made much smaller in magnitude. 

Instead, we propose an alternative form of softening, 
still based on the 2D kernel softening discussed above, with 
the modification that the effective mass m is multiplied by 
an additional modification factor so that 



(3u-|w 2 )m ifu<2/3; 
m if v > 2/3. 



(21) 



The choice of this exact form of softening is arbitrary, but is 
motivated by three desirable conditions on the force within 
and outside the softened region. We require that the force 
decreases linearly to zero at zero separation, yielding a truly 
collisionless form. Second, at separations > 2h, we require 
the force returns to its correctly normalized, perfectly New- 
tonian form and, finally, the derivative of the force is con- 
tinuous, so that the force varies smoothly at all separations. 
Coincidentally, for v < 2/3, this form duplicates the alge- 
braic form of the derivative of the standard kernel, a fact 
that will prove advantageous for obtaining ratios between 
gravitational and pressure forces near unity in section |2~7I 

Throughout this paper we implement gravitational soft- 
ening using equation ll9l with the variable smoothing length, 
h, as its length scale. We will perform separate series' of 
simulations employing either the original mass distribution 
defined by eauation lliJI and incorporating the kernel in equa- 
tion E| or the modified mass distribution of eauation l21l 



2.6 The relative merits of fixed and variable 
softening 

An important advantage of variable gravitational softening 
is that it allows a modeler to soften the gravitational forces 
on the same length scale used for generating the hydro- 
dynamic quantities. The variation is required because the 
hydrodynamic quantities are generated using an approxi- 
mately (or exactly) fixed number of 'neighbors' but the lo- 
cal particle density is not constant. In order to retain the 
approximately fixed neighbor count, the smoothing must be 
correspondingly smaller in high density regions than in low 
density regions. In order to retain the equality over the du- 
ration of a simulation, the softening must also be allowed to 
vary according to local conditions. 

The advantage of softening with the same length scale 
as smoothing is that large imbalances in the gravitational 
and hydrodynamic forces between pairs of particles cannot 
develop due to particles being in range of one or the other 
of the lengths, but not both. The disadvantage is that en- 
ergy may no longer be conserved because no account is made 
of the change in the internal mass distribution of particles 
as their smoothing length changes. Contributions to grav- 
itational potential energy dependent on such changes will 
therefore not be evaluated correctly. 

In the case of fixed softening, an advantage is that it 
conserves energy, with disadvantages including the fact that 
the smallest resolvable length scale is both fixed at the be- 
ginning of the simulation and is the same everywhere. A col- 
lapsing or expanding body may quickly reach a size where 
the flow is dominated by the softening in one region, while 
in another the flow may become unphysically point-like. 

Several previous works have discussed the consequences 
of the alternatives in simulations of different types, but no 
clear consensus has emerged from the discussion. For exam- 
ple, BB97 note that the Jeans wavelength, cast in the form 
of a Jeans mass, must be well resolved by both the gravita- 
tional softening length and by the SPH smoothing length, 
used respectively to ensure numerical stability in the code 
and to produce hydrodynamic quantities from the particle 
distribution. For simulations in which the Jeans wavelength 
is much larger than either the smoothing or softening lengths 
(i.e., that the region is stable against gravitational collapse), 
BB97 claim that little difference in behavior should be ex- 
pected. However, for a marginally stable mass distribution 
(e.g. one Jeans mass distributed over a volume of radius one 
Jeans length) , fragmentation could be artificially suppressed 
or enhanced by changing the ratio of the gravitational soft- 
ening to SPH smoothing length, because of the large force 
imbalances that develop. To avoid such artificial results, they 
recommend a variable gravitational softening set to the same 
length scale as the hydrodynamic smoothing. In the alterna- 
tive, when fixed softening must be used, they recommend a 
softening length no smaller than the value that would cause 
the resolution condition of equation 0] to be violated. 

Other work l|Thacker et al. I l200Cf) fixes the softening 
length, but allows the smoothing to vary, down to a limit 
of either e/2 or e/20. They show that a simulation includ- 
ing cooling will produce much more fragmentation in the 
latter case, and conclude that smoothing lengths must be 
restricted to be greater than the softening length to avoid ar- 
tificial fragmentation. Although implementing fixed soften- 
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ing rather than variable, their conclusion appears consistent 
with that of BB97, but takes no account o f artificially sup- 
presse d fragmentation. On the other hand. Iwilliams et al. I 
( 2004) examine fixed softening models where the smoothing 
is allowed to vary with and without the constraint that it 
may not decrease below the softening length. They conclude 
that the smoothing does not need to be equal to softening 
and that smoothing must not be constrained because hydro- 
dynamic shocks cannot be properly resolved on size scales 
smaller than the softening length. 

2.7 The limits and meaning of resolution in the 
context of softening and smoothing 

For SPH, where hydrodynamic quantities are derived from 
interpolations between pairs of particles, a necessary math- 
ematical condition for the interpolation kernel is that it be 
continuous and have a continuous first derivative (see e.g. 
iMonaehanl Il992l . for additional details). Physically, the re- 
quirement is equivalent to the statement that contributions 
to hydrodynamic quantities like mass density and mutual 
pressure forces undergo no discontinuous jumps as pairs of 
particles approach each other. An unfortunate consequence 
of the two continuity requirements is that for kernels like 
equation 1171 particles that approach each other experience 
pressure forces that first increase with decreasing separa- 
tion, then decrease to zero as the separation between them 
decreases to zero. 

The consequence is unfortunate first, because, to the ex- 
tent that we can regard the approach of two particles (rather 
than a large sample of particles) as representing a compres- 
sion, the fact that the pressure force decreases to zero as they 
approach coincidence means that our intuitive expectation 
that pressure forces continue to increase during a compres- 
sion is violated. Secondly, and as noted above, it may l ead to 
unphysical particle clumping. As described bv lHerant] jl994h 
for non-self gravitating simulations, the reason is that during 
the natural course of a simulation, particles that find them- 
selves closer than a critical separation distance (v = 2/3 
for the kernel in equation II 711 where the pressure force is 
at its maximum, experience relatively smaller mutual pres- 
sure forces pushing them apart, and relatively similar exter- 
nal forces perturbing their motion. Because of these small 
forces, particles continue to travel together with an end state 
in which pairs of particles actually coincide, defining what 
Herant calls a 'pairing instability'. 

The astute reader will not fail to note that discussing a 
'pairing instability' in the context of a SPH simulation would 
seem to be rather irrelevant. SPH particles are of course not 
assumed to represent actual physical particles at all, but 
rather some (statistical) realization of an underlying physi- 
cal system. Therefore any such discussion would appear to 
be meaningless. In our defense, we point out that our discus- 
sion is of a failure mode for the method where the statistical 
assumptions break down, and therefore will require a careful 
analysis. 

As in the case for smoothing, softened mutual gravita- 
tional forces first increase in magnitude as pairs of particles 
approach each other, then decrease to zero as they approach 
still further. In this case however, the decrease is not due to 
any constraint on the kernel, but rather on the combined as- 
sumptions that the softening represents a mass distribution 



and that, at some separation, the mass enclosed by a sphere 
whose radius equals that separation is less than the total 
mass. Then, by invoking Gauss's law, only the fraction of 
mass inside the sphere contributes to the force. In 2D simula- 
tions, where the geometric arguments leading to the Gauss's 
law r esult do not hold (see e.g. figure 3 of iNelson et al. I 
1998), it is necessary to relax the strict identification of the 
kernel with a mass distribution. For the purpose of avoiding 
infinite force contributions at coincidence, softening accord- 
ing to this procedure is effective, though as we show below, 
problems remain due to the fact that the force still does not 
decrease to zero at coincidence. 

In both cases, the reason for the decrease in force magni- 
tude at small distances is to avoid numerical pathologies: for 
gravity, the development of unphysical point-like force con- 
tributions, for hydrodynamics, unphysical discontinuities in 
the forces and other hydrodynamic quantities. The impor- 
tant point to note regarding both softening and smoothing is 
that on such scales, the magnitudes of the forces are depen- 
dent on assumptions made outside the realm of the physical 
model. In other words, the forces computed there are under 
resolved and any developing phenomena sensitive to effects 
in that region can be due only to an external assumption 
rather than to any physical process. 

In order to ensure physically valid simulations, it will 
therefore be important to ensure that effects originating on 
unresolved scales do not drive the results. Moreover, be- 
cause they act on similar scales but with effects of opposing 
sign, it will be important to consider the resolution limits 
of both softening and smoothing lengths together. Setting a 
softening length much smaller than the smoothing length is 
equivalent to the statement that at some distances pressure 
forces are under resolved (i.e. that an error is made in evalu- 
ating them), but gravitational forces are not under resolved 
(that no error is made in evaluating them). The same state- 
ment is true in reverse when the smoothing length is smaller 
than the softening. 

Regardless of whether or not one or the other force actu- 
ally can be resolved better than the other, the numerical as- 
sumptions always limit the effective, physically correct force 
resolution to the larger of the two scales. As pointed out by 
BB97, the consequences of differing values of softening and 
smoothing lengths are net force imbalances of up to a fac- 
tor seven at small separations, even when the softening and 
smoothing are only different by a factor of two. Much larger 
force imbalances develop when the ratios become more ex- 
treme, and BB97 attribute the unphysical outcomes of sev- 
eral simulations presented in the literature to this source. 

Figure shows the ratio of the gravitational and pres- 
sure forces between two particles located in a disk where 
the Toomre Q value set to unity, realized in 2D. As for 
the 3D case, forces are nearly balanced when softening and 
smoothing are set to identical scales, but large imbalances 
are present when they are not identical. Unlike the 3D case 
however, the force ratios for the standard kernel soften- 
ing/smoothing option do not remain near unity as the sep- 
aration, v, goes to zero, but instead approach infinity. Due 
to the reduced dimensionality, the gravitational force ap- 
proaches a finite, non-zero value rather than decreasing to 
zero as it does in the 3D case. Force imbalances such as these, 
whether caused by unequal softening and smoothing lengths, 
or by the form of the softening or smoothing themselves (or 
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Figure 1. Ratios of the gravitational to the pressure forces for 
the three scenarios we test, each as a function of the normal- 
ized separation, v = r/h, between two particles. The top panel 
shows ratios using the standard kernel derivative for the pressure 
force (top) and the standard kernel softening for the gravity. The 
middle panel shows the ratios using pressures obtained with the 
alternate kernel derivative of equation ll8l with standard softening. 
The bottom panel shows ratios using the modified kernel soften- 
ing with the standard kernel derivative. The solid lines in each 
panel refer to the condition where the softening and smoothing 
lengths are equal, while the dotted and dashed lines denote the 
condition where the softening is twice and half of the smoothing, 
respectively. 



both), may influence the susceptibility of the simulation to 
artificially induced fragmentation. When the imbalance fa- 
vors gravity, particles are not only passively drawn to each 
other, but are also actively attracted. Given an imbalance 
of the opposite sense, particles may be artificially repelled 
from each other, preventing fragmentation of physical origin 
from occurring. 

Given the intrinsic force imbalances present in 2D sim- 
ulations using the standard kernel for smoothing, even 
when used with identical length scales, simulations may 
therefore be susceptible to artificially induced fragmenta- 
tion. On the other ha n d, the modified kernel derivative of 
iThomas fc Couchmanl (^992) allows the force ratio to re- 
main near unity all the way to zero separation because the 
pressure approaches a finite, non-zero value in step with 
softened gravity. Similarly, the modified gravitational soft- 
ening defined by equation also results in force ratios near 
unity for equal softening and smoothing lengths, but in this 



case, both forces approach zero as their limiting value rather 
than a finite quantity. In both cases, large force imbalances 
can develop when the softening and smoothing are unequal. 
With the modified softening variant, force imbalances of up 
to a factor ~ 7.5 at v = develop when the lengths are 
not equal, nearly twice the imbalance present with the al- 
ternate kernel derivative. Because the actual magnitudes are 
smaller, the net influence of the imbalance might be smaller. 

We will investigate this question, the influence conse- 
quences of using the alternate kernel derivative formulation 
and the alternate softening formulation and the importance 
of using fixed softening or variable softening for which force 
imbalances may or may not develop respectively, in our work 
below. 

2.8 The relative merits and shortcomings of 2D 
and 3D treatments of the disk 

As computational facilities have become more and more 
powerful, simulations of greater and greater complexity have 
become possible to perform at acceptable cost. An important 
step in increasing complexity is to increase the dimensional- 
ity, first from one to two dimensions, and then to three. In 
the context of circumstellar disks, the full transition to 3D 
remains incomplete. Some workers prefer to limit their sim- 
ulations to two dimensions, while others have begun work in 
3D. A number of significant consequences derive from each 
choice. 

Of primary importance is the computational cost and 
its relation to the spatial resolution affordable to the simu- 
lation. For a grid based method, computational cost in 3D 
will be approximately proportional to the fourth power of 
the number of cells in any one dimension, while for 2D, the 
proportionality drops to the third power. A similar propor- 
tionality will hold for particle simulations as well. Therefore, 
simulations in 3D will be able to employ fewer total cells 
or particles per dimension: linear resolution is intrinsically 
coarser in 3D than in 2D for simulations of comparable cost. 
In contrast, while 2D simulations allow dramatically higher 
linear resolution in the two dimensions they actually model, 
they do so at the cost of requiring assumptions be made 
about the character of the system and its behavior in the 
third dimension. 

For a 3D simulation, if it is to be truly three dimen- 
sional, the physical extent the system in each direction must 
be resolved by some number of particles or grid cells greater 
than one. The exact minimum number will be a function 
of both the problem and of the method employed to evolve 
the system. In the context of circumstellar disks, the fact 
that the disk is spatially thin means that the available res- 
olution must be allocated inhomogeneously in space if the 
cost of calculation is not to become too exorbitant. In a grid 
simulation for example, many more grid cells must be allo- 
cated to the same physical length in the z coordinate than 
either the x or y coordinates, to resolve the vertical extent of 
the disk. Problems may arise from asymmet ric grid resolu- 
tion, including incorrect gravitational forces (iPickett et al. I 
2003), or incorrect hydrodynamic evolution. 

The problem is especially acute for SPH simulations 
because the kernel used to reconstruct the hydrodynamic 
quantities is spherically symmetric. Although experiments 
with non-spherical kernels have been published, they are 
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not common, due in part to the difficulty of implementing 
them in a form that retains angular mo mentum conserva- 
tion fe.g. lFulbright. Benz fc Davieslll99BF) . With a spherical 
kernel and in the low resolution limit, particle smoothing 
lengths may actually exceed that of the disk's vertical ex- 
tent, so that only a single fluid element spans the entire 
system in that coordinate. Similar conditions may occur in 
grid based simulations, if extremely asymmetric zone dimen- 
sions are not to be encountered. Any 3D simulation carried 
out under such conditions will effectively model only two di- 
mensions, but will not include any supporting assumptions 
present in a simulation explicitly limited to 2D. 

One such assumption will be in the description of the 
fluid itself. In most 'normal' fluids, pressure is by definition 
a scalar quantity whose gradient causes a force to be exerted 
in all directions. The procedure used in SPH to construct the 
hydrodynamic quantities however assumes a roughly spher- 
ical distribution of particles. Otherwise, the interpolations 
at the heart of SPH are no longer interpolations but in- 
stead extrapolations in directions where few particles exist. 
As a result, and as is well known to its practitioners, the 
method becomes very inaccurate at boundaries. In circum- 
stellar disk simulations, performed with SPH at resolution 
where only one or few particles span the entire vertical ex- 
tent of the disk, essentially all particles will be located near 
boundaries. 

Whether or not a particular simulation with one dimen- 
sionality will be more physically meaningful than a compara- 
ble one of the other depends on the quality of the assump- 
tions made about the third dimension in a 2D simulation 
and whether those assumptions offset the loss of linear res- 
olution affordable in 3D. In sectional we will show that 3D 
simulations of disks require much higher resolution than one 
might naively expect in order to reproduce basic hydrody- 
namic features of the flow correctly. Such high cost means 
that a given study will be able to perform many fewer sim- 
ulations, severely constraining the possibility of performing 
the large parameter studies often required to fully explore 
the implications of a given physical model. We have there- 
fore concentrated on the study of disks in 2D throughout 
the rest of this paper. 

2.9 The interpretation of hydrodynamic 
quantities and gravity in 2D 

Because disks are in fact truly three dimensional, in spite 
of our approximation that they are thin, it will be impor- 
tant for correctly interpreting the results of the simulations 
to understand the consequences of a 2D approximation and 
where it may break down. For 2D simulations, two funda- 
mentally different assumptions about the third coordinate 
are possible. The modeler may either assume that the simu- 
lation is modeling an infinite cylinder or that the simulation 
is modeling a thin system in which the system's dynamics 
and morphology are either negligible in the third dimension 
or some approximation is made regarding their behavior. 
Each assumption leads to quite different treatments of the 
hydrodynamics and gravitation in the simulation. 

In the case of the infinite cylinder assumption, each 
point in the plane actually corresponds to a line extending 
to infinite distance in the positive and negative third coordi- 
nate, which we will assume to be the z coordinate. For sim- 



ple hydrodynamical problems, the interpretation poses no 
particular conceptual difficulty since hydrodynamic quanti- 
ties will carry over from 3D to 2D unaltered. The conse- 
quence for gravitational or electrostatic forces however, is 
that they become inversely proportional to the separation 
in the xy plane, rather than to the inverse square of the 
separation. For a thin system, the opposite situation holds. 
Inverse square law forces retain their familiar 3D form, but 
hydrodynamic quantities must be altered, specifically into 
integrals of their true 3D forms. For example, surface den- 
sity may replace volume density. 

The circumstellar disks in this study are spatially thin, 
and the most natural interpretation of a 2D simulation of 
such a disk is that of a thin, vertically integrated model. 
In keeping with this characteristic and with many previ- 
ous analytic and numerical treatments of disks, our 2D disk 
simulations are performed in this context. In order to allow 
readers to evaluate the results of our study more thoroughly, 
we now discuss several factors important for their interpre- 
tation and similar work by others. 

In addition to the requirement for 2D simulations that 
state variables be integrated over the z coordinate, a more 
subtle modificatio n must also be made to other hydrody - 
na mic quantitie s. |GoJdreidL_G podman fc Naravanl (ll986T> 
and lOstriker. Shu fc ATkamsTll992l) each discuss the mod- 
ifications in the effective value of its polytropic index, n, of 
the gas, corresponding to a degree of freedom correspond- 
ing to the disk 'puffing up' in the third coordinate. Both 
conclude that a value of 7 slightly reduced from that ex- 
pected for the 3D case should be used, due to the additional 
freedom. When an isothermal equation of state is employed, 
the 7 value will retain its limiting value of unity, so no affect 
will be present in the work here. The change will also not 
be required in a truly 'razor thin' 2D model, where motion 
is restricted in the third dimension entirely. In this case, the 
effective 7 would increase instead, due to a decrease in the 
number of internal degrees of freedom for the gas. 

Vertically integrated quantities also require special con- 
sideration in realization and interpretation of gravitational 
forces in the system, even though they retain their familiar 
inverse square form. As noted in section l2.5.2l gravitational 
forces obtained from a straightforward carryover of the 3D 
method of kernel softening exhibit a non-zero force at zero 
separation. From a physical perspective, such a condition 
will simply be erroneous, since in reality the disk mass de- 
scribed by the particles is spread over some vertical extent: 
it is no longer 'thin' compared to interparticle separations. 
Where in reality, the force of one vertical column on an- 
other falls to zero at zero separation, the force based on a 
vertically integrated mass located at the disk midplane does 
not. 

From a numerical perspective, this condition challenges 
the assumption that the particles are collisionless, in turn 
the assumption that softening was introduced to ensure. 
While considerably weakened, consequences will be less se- 
vere then in the 3D case because the interaction force does 
not become infinite at any separation, as it would if two 
truly collisional particles were to interact. We may therefore 
attempt to salvage the gravitational forces via some simple 
modification, as we suggest with the altered softening pre- 
scription defined in equation 1211 or the altered description 
of pressure forces derived from equation 1181 
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Care must still be taken using either modification in any 
2D simulation. Interparticle forces will deviate from their 
vertically integrated forms over spatial scales comparable 
to the scale height of the disk. At extremely high resolu- 
tion, where particles correspond to very thin columns of fi- 
nite extent in the vertical direction, the softening or ker- 
nel derivative modifications will affect the calculated forces 
only over a small fraction of that distance. At separations 
between the particle size, h, and the scale height, H, forces 
will therefore be overestimated. Numerical experimentation 
has shown that deviations will not be large at resolutions 
where h>0.1H, but become much more significant when 
h ~ 0.01//. Simulations presented here do not fall into the 
latter category. 



3 TEST PROBLEMS FOR DETERMINING 
THE REQUIRED RESOLUTION OF 
SIMULATIONS OBEYING THE TOOMRE 
CRITERION 

What resolution is required (i.e. what values of T and N leBO 
from equations 151 and I1H to ensure that a simulation that 
produces collapsed objects in a disk is producing numeri- 
cally valid results? As was done to develop criteria for Jeans 
collapse, we will define a specific problem on which to com- 
pare the results of a numerical code at different resolutions. 
In parallel, we will investigate the influence of three different 
strategies for the treatment of small scale interactions be- 
tween particles: a 'base' version and two variants that mod- 
ify the kernel derivative in one case or the kernel softening 
in the other. 

We use a small variation ( see below) of a simulation 
discussed in lNelson et al. I lll99ct) . who used an SPH code to 
model the evolution of disks in two dimensions. Simulations 
using SPH are especially sensitive to violation of a resolu- 
tion criterion because resolution is dynamically allocated. 
Particle smoothing lengths are ordinarily considered to be 
functions of the local flow variables so that in high density 
regions, they shrink in an attempt to follow the small scale 
motions of the fluid there. In most respects, this feature can 
be extremely desirable because there is no a priori reason 
to expect fragmentation in one or another part of a given 
simulation. On the other hand and as we show below, in- 
sufficient care in its use can lead to numerically induced 
fragmentation. 

We use the VINE code (Wetzstein et al. , Nelson et 
al. , in preparation) in its 'SPH only' mode and using its 
leapfrog integrator to perform our simulations. An earlier 
version of this code, with a second order Run ge-Kutta inte- 
grator, was used in the original calculations in lNelson et al. I 
1998. Exploratory tests with VINE using this same integra- 
tor in the present simulations showed that similar results 
were obtained from both. Most calculations were performed 
with a single, global time step for all particles in order to 
assure that our results were unaffected by as few systematic 
effects as possible. Some tests were made with individual 
time steps in order to explore effects of numerical stabil- 
ity due to this source. These calculations required ~ 3 — 5 
times less computer time to complete and, while some differ- 
ences between global and individual time step versions were 
present in the results, none materially affect the conclusions 



made from them. VINE uses a binary tree to organize par- 
ticle data, so that they may be accessed efficiently for use 
in both the hydrodynamic and the gravitational force calcu- 
lations. In order to avoid calculation times for the gravita- 
tional forces of order 0(N 2 ), sufficiently distant particles are 
approximated as nodes in the tree, resolved to quadrupole 
order in the actual calculation. The acceptability criterion 
for the nodes was set so that forces on > 99% of particles 
would be accurate to <0.1%. 

VINE employs an artificial viscosity with both bulk 
and von Neumann-Richtmyer terms to stabilize the evolu- 
tion and to convert kinetic energy into thermal energy in 
shocks. The coefficient for each term were set to a = 1 and 
= 2, which are values standard in the literature. Using 
these values, simulations of disks using SPH are afflicted 
with a quite large and unphysical level of shear viscosity. In 
order to minimize such effects, the simulatio ns here were ru n 
with the shear viscosity reduction switch of iBalsaral il995t) . 
which reduces the magnitude by a factor ~ 3 — 5. 

Important parameters from the simulations presented 
here and in the following sections are listed in Table [3f. The 
columns of the table show the resolution of each simulation, 
the type of gravitational softening used (for SPH simula- 
tions) and, finally, the time at which the time at which the 
first clump is formed and the duration of the simulation. 
Simulation names ending in 'TC and 'grv' denote simula- 
tions run with the modified kernel derivative of equation 
1181 or the modified gravitational softening of equation 1211 
respectively, while those suffixes represent identical initial 
conditions but run with the unaltered kernel derivative or 
softening. All simulations are performed in 2d and include 
the effects of self gravity, except the sgoff2d4 simulation in 
2D which does not, and the remaining sgoff simulations, dis- 
cussed in section^] which also do not include self gravity and 
are done in 3D. 



3.1 The definition of our test problem 

iNelson et al. I (jj.998) modeled the evolution of self gravitat- 
ing disks in 2D with masses between 0.05 and 1.0 times the 
mass of the central star, using SPH. Other simulations from 
that work, performed using PPM, are not considered here 
because they could not be carried out far enough into the 
high amplitude regime to make fragmentation likely. The 
specific model we consider here (labeled 'scv2' in that work), 
had an assumed disk mass of Md=0.2M 4 and a minimum 
Toomre Q value defining its stability of Qmin =1.5. We be- 
lieve that model will be a particularly challenging test of the 
criterion because collapse was observed after only about one 
orbit of the outer disk edge, corresponding to about 11 or- 
bits in the region where clumps first started to form. By the 
conclusion of the run at 1.6Tb, more than 30 clumps had 
formed. As in the originals, we have simulated the evolution 
in two dimensions, so that the surface density is directly 
available from the calculation. 



2 Note to MNRAS latex programmer/copy editor: It would be 
nice if the reference to the table automatically came out correctly: 
There is only one table in this paper, but latex calls it table 3 
here and elsewhere in the text, and table 1 in its definition. 
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Table 1. Simulation Parameters 



Label 


Resolution 


Softening 


T first 


Tend 


modi 


7Q/I A 


var. 


U.ooi /) 


1 9 Tr, 
l.Z i/) 


ITIOCIZ 




Var 

v ar. 


1 n7T^ 


1 fi TV 
1.0 I/) 


mod3 




Var 


^ Q1 


D.U 1 d 


inod4 


9^091 
ZOUZlo 


Var 

v ar. 




in n t~" 

1Z.U 1 £) 


mr»/-] 1 TP 


7QzlA 


Var 

v ar. 


1 1 

1.101/j) 


i Tr, 

1.0 lj) 


moaz i 


Q91 99 


Val- 
var. 


1 OfiT„ 

l.zol ]j 


1 7^Tr, 
1. f Oi_D 


moQo i 




Var 

v ar. 




in n a - * 

1Z.U 1 


mr.fl/ITP 

lnOQ4 1 \^ 


9R091 Q 


var. 




1 9 n t», 

1Z.U 1 £) 


fivl TP 


Q91 99 


nx.x. A TT 


u.uoi x> 


n 1 7Tn 

U.l ( 1 £) 


fiv9TP 

nxz i o 


QOI 99 
OZ±ZZ 


Afl ATT 


9 1 ^TV. 
Z. 101 JJ 


n Tr, 

O.U 1/9 


fix3TC 


32122 


fi 1 ^ ATT 




^ n Tr, 

O.U 1 £} 




Q91 99 


11 ATT 




c n t 

O.U Ji9 


m 1 fVTP 

rn i ix i kj 




09^ ATT 
.UZO-rVU 


o ^oTt. 

U.OZl £j 


U.O I 1 £) 




ZOUZ lo 


14 ATT 
. 14 AU 




19 Tr, 

1Z.U 1 £) 


modlgrv 


7QAA 


v ar. 


u.yoi jj 


1 97Tr, 


mod2grv 


Q91 99 


val- 
var. 


9 OST^ 
Z.Zol /J 


9 Tr, 
Z\0 ±D 


mod3grv 


1 9QQQ/1 


Var 




19 Tr, 
1Z.U 1 £) 


mod4grv 


9fifl91 
ZOU/lo 


Var 

v ar. 




19 Tn 

1Z.U 1 £) 


fixlgrv 


QQl 99 


o^ ATT 


U.U i 1 JJ 


90T,-, 
U.ZU1 £) 


fix2grv 


QQl 99 


AO ATT 


1 A9Tt^ 


1 Q Tr, 

i.y l ]j 


fix3grv 


32122 


fi^ ATT 




^ Tr, 

O.U 1 £) 


fix4grv 


32122 


1 1 AU 




^ Tn 

O.U 1 £) 


rn 1 fxgr v 


9R091 ^ 
ZOUZlo 


09^ A TT 


fl 9QTr, 

u.zyi £j 


o f;nT n 


m4fxgrv 


260213 


.14 AU 


2.20T D 


2.50Td 


sgofT2d4 


260213 


Var. 




12.0 T c 


sgoff3 


129384 


Var. 




2.0 T D 


sgoff4 


260213 


Var. 




2.0 T D 


sgoff5 


502089 


Var. 




2.0 T D 


sgoff6 


994740 


Var. 




2.0 T c 


Boss 


100x23x256 


Grid 


345yr 


359 yr 



The temperature and surface density of the gas were 
defined in that model with softened power laws as: 



T(r) = To 



and 



E(r) = So 



1 + 



1 + 



-q/2 



-p/2 



(22) 



(23) 



where q = 1/2 and p — 3/2, respectively, and the core radius 
for both power laws was set to r c = 1. The constants To and 
Eo were determined from the assumed Toomre stability pa- 
rameter and the radial dimensions of the disk, defined at its 
inner edge by Ri — 0.5 and its outer edge by Rd = 50. Par- 
ticles were laid out on concentric rings and an equilibrium 
state accounting for stellar and disk self gravity as well as 
pressure forces defined the velocities of each particle. The gas 
was evolved under the forces of stellar gravity, self gravity 
and gas pressure, which was computed using an isothermal 
equation of state (i.e. 7=1 with fixed temperature as a 
function of radius) . Because of the simple physical assump- 
tions, the dimensions of the system were scalable. Given a 
star of one solar mass and a disk radius of 50 AU, one orbit 
at the outer edge of the disk requires approximately 353 yr 3 . 



3 Note that these u nits are slightly different than those used in 
iNelson et al. lfl998fl . where a star of 0.5Mq was assumed in order 
to correspond more closely to typical T Tauri star masses. 



During the preparation of this manuscript, we deter- 
mined that simu lations evolved using initial conditions iden- 
tical to those in lNelson et al. I il99ct) were not possible be- 
cause at high resolution the large pressure gradient at the 
disk's outer edge caused unphysical behavior in the system. 
In order to sidestep this problem, we modified the form of 
the surface density power law near the outer disk edge so 
that the discontinuity is spread over a larger radial range. 
In this work, the surface density is distributed according to: 



E mod (r) = SE(r) 



(24) 



where the factor, S, is a linearly decreasing function near 
the disk boundary and is defined by: 

(I for r < R D - 8 ; 

S = < 1 — (rd&2=m for R D - 8 <r<R D +8; (25) 
{ " for r > R D + 8 ; 

With this definition, the disk edge is smoothed in the region 
within a distance 8 inward and outward of the nominal disk 
radius. In the simulations here, we define the smoothing pa- 
rameter 8 — 5 AU. The other initial conditions and physical 
assumptions remained the same. 

3.2 Evolution of our test simulations 

Using the conditions defined above, we ran a set of four sim- 
ulations with varying resolution using each of three treat- 
ments for the small scale interactions between particles, but 
otherwise identical. The first treatment, with simulations la- 
beled modi-mod^ in tabled] uses gravitational softening of 
the form defined by equation 1191 and the standard kernel 
gradient defined derived directly from equation ll7l The sec- 
ond treatment, with simulations labeled modlTC-rnod^TC, 
uses the modified kernel derivative defined by equation 1181 
to determine mutual pressure forces between particles, again 
with the standard kernel softening. The third treatment, 
with simulations labeled modlgrv-mod^grv, uses the modi- 
fied gravitational softening defined in equation 1211 with the 
standard kernel gradient. Figures |3 El and 0] show these four 
realizations of the model, each at the termination of the 
simulation. 

Each model develops multiple armed spiral struc- 
tures over most of its radial extent, as in the previous 
INelson et al. I |l998) work. The spiral structures are filamen- 
tary and change the details of their appearance as the sim- 
ulation proceeds. Evolution after the formation of the spiral 
structure however, was strongly dependent on the resolution 
employed. In evolution of the three lower resolution real- 
izations of the mod series, spiral struc tures eventually frag- 
mented into multiple clumps, as in the INelson et al. I (1998) 
work. The time at which clumps begin to form, measured 
from the beginning of the simulation, was later at higher 
resolution than at lower. Moreover, fewer clumps formed in 
the higher resolution realizations than in the lower, and in 
the highest resolution realization, no clumps were produced 
at all: the delay before clumps begin to form has become 
longer than the duration of the simulation (12Td, or about 
4200 yr). 

Similar statements are true of both the modTC and 
modgrv runs, though there are a number of important dif- 
ferences as well. Of particular interest is the fact that the 
structure seen in figures [3] and 0] is substantially smoother 
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Figure 2. The particle distribution for the mod series of simulations at the end of each simulation. At the times shown, filamentary spiral 
structures have developed throughout the disks and are visible as azimuth varying blue, green and yellow color variations in the images. 
The three lowest resolution realizations have also produced a number of clumps (visible as red dots in the image) typically containing 
as many as several hundred individual SPH particles. Clumps in the mod3 realization are visible at about the 9 o'clock position in the 
figure, while in the two lowest resolutions, they are distributed throughout. The highest resolution realization produced no clumps. The 
color scale shown defines the base 10 logarithm of the surface density. 



than in the corresponding panels of figure[5] The overall rel- 
ative smoothness is also reflected in the number of clumps 
that form in each realization. Where the modS simulation 
formed several clumps, simulations with the two modified 
treatments modSTC/modSgrv did not. Although each of the 
corresponding lower resolution realizations did form clumps, 
the number that formed was smaller in each case and delayed 
in time relative to the standard kernel versions. 

Figures Q>] |S| and |7| show azimuth averaged, radially 
binned surface density profiles of the same models as shown 
in figures |2] EH and 01 at the same times. Also shown are 
instantaneous maximum surface densities, defined as the 
maximum in each ring seen at the time shown. The aver- 
age in each bin is weighted by the number of particles in 
the bin, with the width of each bin set to 0.02 AU. At small 
radii, accretion onto the star has depleted the initial power 



law density distribution, in each case by a differing amount 
due to the differing magnitude of the dissipation derived 
from the artif i cial v iscosity included in the simulation (see 
iNelson et al. 1 12000) for a discussion of this effect). Further 
out, a number of spikes are visible in the profiles, each cor- 
responding to the radial location of a clump in the disk. The 
lowest resolution realizations display a very large number of 
clumps, while successively higher resolution realizations pro- 
duce fewer, or none. Both the number of clumps that formed 
and the radial extent over which they are found are smaller 
when we use the two modified kernel treatments, compared 
to the standard form. Though not possible to see in snap 
shots of a single time, an important difference between the 
evolution with the standard and either of the modified ker- 
nel treatments is that the cumulative maxima in the latter 
case remain far lower than the former and also do not appear 
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Figure 3. The particle distribution for the mod series of simulations with the TC92 kernel derivative, each at the end of each simulation. 
As in figure [5] filamentary spiral structures have developed throughout the disks and are visible as azimuth varying blue, green and 
yellow color variations in the images, with clumps (in the two lowest resolution realizations) visible as red dots in the image. The color 
scale is the same as in 121 



to increase over time. We believe that this near steady state 
behavior represents the correct evolution of the system, even 
if it were to be evolved further in time. 

What is the origin of the differences in the morphol- 
ogy between each series and the others? We postulate that 
it lies in the treatment of the small scale interactions be- 
tween particles, and specifically in force imbalances between 
pressure and gravity that may be present. If so, we may ex- 
pect distortions in the distribution of interparticle separa- 
tions, relative to what may be expected and to what may be 
present with other realizations. Figure [S] shows histograms 
of the distances between all neighbor particles for each of 
the three realizations of the mod4 simulations, as well as a 
simulation with the same resolution and initial condition, 
but without self gravity. Each histogram bin is defined to 
be of width 5v — 0.002, and all neighbors for all particles 
in each simulation are accounted for in the histograms. On 



the scale of the smoothing length, we can approximate the 
surface density as constant. For a uniform density distribu- 
tion in a 2D simulation, realized with a random distribution 
of particles, we expect that the number of neighbors to be 
found at a given separation to increase linearly with the 
separation itself. At large separations, we recover the linear 
behavior. At smaller separations, all four distributions de- 
part from the linear proportionality because the assumption 
underlying the linear proportionality neglects the influence 
of interparticle forces. Due to the different treatments of the 
small scale interactions, some differences also appear in the 
neighbor distributions. 

The most striking feature of the plot however is the 
behavior near v = 0. In the standard large num- 

ber of particles actually coincide in space, while in each 
of the other three variants, none do. This coincidence rep- 
resents an extreme example of Herant's pairing instability 
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Figure 4. The particle distribution for the mod series of simulations with the alternate gravitational softening, each at the end of each 
simulation. As in figure [2] filamentary spiral structures have developed throughout the disks and are visible as azimuth varying blue, 
green and yellow color variations in the images, with clumps (in the two lowest resolution realizations) visible as red dots in the image. 
The color scale is the same as in 121 



discussed above, that develops due to the interparticle pres- 
sure vs. gravitational force imbalance at small separations. 
The pairing instability seen here is clearly not identical to 
that discussed by Herant however, since no pairing is present 
in the non self gravitating realization for which conditions 
are most similar to those discussed by Herant. The pair- 
ing is also not similar to t he numerical instabilities seen by 
llmaeda fc Inutsukal i2002tl for the same reason. 

Figure[§]shows that the number of paired particles con- 
tinues to grow as the simulation proceeds. After 12Tb, more 
than 65000 pairs have formed from > 130000 particles: of the 
~ 260000 original particles, about half have become paired. 
The existence of paired particles means first that the ef- 
fective resolution of the simulation decreases with time as 
more and more particles become paired, and is consistent 
with our qualitative observation above that the cumulative 
maxima in these simulations also increased. Their presence 
also means that even though the force contribution from a 
single neighbor particle is small relative to the individual 



contributions from the rest of the system, when the par- 
ticles approach coincidence the pairwise contribution can 
have an important and large scale influence on the behav- 
ior of the simulation. Based on these results, and the im- 
pact they have on the large scale disk morphology and its 
tendency to fragment, we conclude that using the standard 
kernel based gravitational softening in combination with the 
standard kernel derivative is likely to produce results con- 
taminated with numerical artifacts in 2D simulations, and 
should not be used. 



3.3 Determining the resolution criterion 

For each of the three series' of simulations, the propen- 
sity of simulations to fragment decreased as resolution in- 
creased, until at sufficiently high resolution no fragmenta- 
tion occurred at all. This behavior reflects that seen by 
iTruelove et al. I (119971) . where fragmentation was enhanced 
when resolution was insufficient, rather than that seen by 
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Figure 5. The azimuth averaged and in stantaneous maximum surface densities plotted as functions of distance from the star for the 
test model derived from a simulation from lNelson et al. I Jl99St) . realized at four resolutions (the 'mod series of simulations— see table 131 . 
each at the end of the simulation. Each panel in this figure corresponds to the same panel as in figure 121 



0x10* 
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Figure 9. The number of particle pairs in simulation mod4 as 
a function of time. Corresponding figure for the other mod4 vari- 
ants are not shown, since no particles become paired during their 
evolution. 



BB97, where fragmentation was delayed. This is fortunate 
because it allows us to use the change in behavior of re- 
alizations of the same initial condition but differing resolu- 
tion to determine empirically the approximate resolution (in 
number of particles) necessary to obtain 'correct' evolution. 
Specifically, we can apply the criterion to two simulations 
which straddle the boundary between those that produced 
fragments and those that did not. We can then note the 
value of iVr OSO at which the criterion succeeds for the entire 
simulation at the higher resolution, but fails for the lower 
resolution version. BB97 scaled the value of the Jeans resolu- 
tion criterion through the value of N Teao as a multiple of the 
average number of neighbors, 7V nc i g h. It will be convenient to 
scale the two dimensional criterion similarly here although 
in both cases, the correct measure will be a quantity inde- 
pendent of the specific neighbor count. We also note that in 
two dimensional SPH simulations, it is usual to use a smaller 
number of neighbors for each particle due to the lower di- 
mensionality. In this work, we have used a number ra nging 
betw een 10 and 30, depending on the local flow (see iBenj 
1990, for details). We therefore scale by a factor iV nc i g h= 20 
for these simulations. 

Figure [TU1 shows the cumulative maximum surface den- 
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Figure 6. The same as figure[F] but for the TC versions of the simulations. Each panel in this figure corresponds to the same panel as 
in figure 131 



sity binned as a function of radius, along with the maximum 
resolvable surface densities determined from equation ^3 us- 
ing three values of N leBO equal to 1, 6, and 12 times the 
average number of neighbors for SPH particles evolved in 
2D (i.e. 20). The cumulative maximum for each given bin 
is defined as the maximum surface density achieved by any 
particle in that bin over the entire course of the simulation 
up until the time shown. 

For the mod3 realization, the cumulative maximum 
exceeds the critical value using the A r rcso = 12iV ne igh crite- 
rion for all radii inside ~ 20 AU. At the same time, the 
A'reso=6A nc igh criterion is not violated except at a few lo- 
calized radii. Interestingly, a small density spike relatively 
early in the simulation near 7 AU did not lead to clump 
formation, but later interactions near 10-15 AU did. For the 
mod4 simulation, the N Teso = 12Ar nc i g h condition is satisfied 
over the entire radial range for the life of the simulation, and 
no clumps formed. We can make only a tentative assignment 
of the value of N leao from this result however because of the 
importance particle pairing may have on the effective reso- 
lution. 

Fieure lTTI shows the cumulative maximum surface densi- 
ties for two of the simulations using the TC92 kernel deriva- 
tive and the modified gravitational softening. In this case we 



plot the curves for mod2TC/mod3TC and mod2grv/mod3grv 
rather than for modS and mod4 in order to maintain the 
straddle of the fragmenting/non- fragmenting outcomes in 
these simulations. In both variants, the cumulative maxi- 
mum densities in the mod2 realizations exceed the critical 
value using the both the N leso =6N n ei s h and A r reS o = 12iV llc igh 
criteria, while the mod3 realizations obey the 6A ne i g h cri- 
terion but violate the 12A\ e igh criterion. Since no clumps 
formed during the evolution, the latter must be considered 
too conservative, and we conclude that the resolution re- 
quired to avoid fragmentation due to unphysical growth of 
self gravitating structures in the disk is A' I eso=6A? n eish- 

At first sight, the required resolution for the Toomre 
condition appears much larger in comparison to that for the 
Jeans condition analysis discussed by BB97. The apparent 
paradox is resolved if we note, as in section T2.2I that their 
assumed value of the Jeans mass was quite low and resulted 
in a relatively lower neighbor requirement. Our definition of 
the Toomre mass corresponds to an analogue of the larger 
definition of the Jeans mass defined as in equation 1131 

The origin of the large required value of iV reS o becomes 
clear when we observe that the values of the cumulative max- 
ima (or indeed, also the slightly lower instantaneous maxima 
not shown) are typically a factor of several higher than the 
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Figure 7. The same as figure[S] but for the alternate gravitational softening versions of the simulations ('grv'). Each panel in this figure 
corresponds to the same panel as in figure HI 



averages. Most of the difference can be accounted for by the 
existence of spiral structure in the disk, however fluctua- 
tions due to the effect of the exact, time varying positions 
of particles relative to each other on the calculation of the 
density make up a smaller, but still significant component 
of the difference. Such fluctuations are intrinsic to the SPH 
method itself and, while a small decrease is observable be- 
tween the low and high resolution simulations in figure [TU1 
their existence will be a part of all SPH simulations. 

It is interesting to note that the densities can exceed 
the resolvable maximum in the lower resolution realizations 
for some time before clump formation begins. Further, the 
amount of time before the onset of clump formation is a res- 
olution dependent quantity, with higher resolution leading 
to more time before clump formation. It is therefore clear 
that resolution studies are particularly important for par- 
ticle simulations showing evidence of clump formation. Ap- 
plying this statement to all three variants of our own mod 
series' of simulations, it is clear that if clumps do indeed 
form in disks similar to those studied, the process requires 
a time scale than longer the 12Tb (~4200 yr) for which we 
have evolved the system. 

In the two lowest resolution realizations of all three vari- 
ants (i.e. modi /TC/ 'grv and mod2/TC/grv), the numerical 



stability criterion is violated even in the initial condition. 
Since these simulat ions are only slight modifications from 
those presented in iNelson et al. I lll998l) and at the same 
resolution, the same statement applies to those simulations 
as well. Regarding the clump formation seen in the simula- 
tions from that work, we originally concluded that although 
the physical model (in particular the isothermal equation of 
state) was insufficient to correctly model clump formation, 
if it was to occur at all, it would be most likely in between 
ab out 10 and 40 AU. Th e conclusion that the physical model 
in lNelson et al. may be insuffi cient may indeed re- 

main valid (see e.g lNelson et al. l2000l) 'l. However, the loca- 
tion and existence of clump formation is definitely invalid: 
it was due purely to failure of the Toomre criterion and not 
to any physical process, however modeled. 



4 FIXED VS. DYNAMICALLY VARIABLE 
GRAVITATIONAL SOFTENING IN 
PARTICLE SIMULATIONS 

After investigating the effects of resolution and choice of ker- 
nel on disk simulations, in this section, we turn to an inves- 
tigation of gravitational softening. We first discuss factors 
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Figure 8. A histogram of the number of SPH neighbor particles found at a given separation, for four variants of the mod4 simulation, as 
labeled. The separation is normalized to the dimensionless variable, v = r/h, and the histograms account for all neighbors of all particles 
in each simulation. 



that motivate the selection of the value of a fixed softening 
length, then perform a side by side comparison of otherwise 
identical simulations employing either fixed or variable soft- 
ening. Due to our conclusion in the last section that the stan- 
dard softening/smoothing combination should not be used, 
we limit our investigations here to the TC92 and modified 
gravitational softening variants. 

4.1 Choosing a fixed softening length 

Unlike the case of variable softening, where the softening 
length is set locally and dynamically, with fixed softening 
the modeler must make a specific choice of the length that 
is appropriate for the entire simulation at all times. Choosing 
an appropriate value in disk models is substantially less triv- 
ial in the case of disk models than in the cloud collapse sce- 
nario discussed by BB97 because, simply as a consequence of 
the initial conditions (i.e. the density gradient as a function 
of distance from the central object), the smoothing lengths 
of particles are relatively steep functions of position. This 
is important because for any constant choice of softening, 
the relative magnitudes of softening to smoothing will also 
vary with position, perhaps artificially suppressing fragmen- 



tation in one region while enhancing it in another. Moreover, 
because the Toomre wavelength will itself be a function of 
radius through both density and temperature, both lengths 
will vary relative to it as well. 

Naively, one might expect that variable softening could 
lead to results that are quite susceptible to small scale frag- 
mentation. For example, if a region begins to collapse and 
the particles representing it approach each other, their mu- 
tual gravitational attraction continues to increase as their 
softening and smoothing lengths decrease, perhaps instigat- 
ing the very problem softening is meant to avoid. It is not 
clear that such a condition exists in practice however. Vari- 
able softening was implemented for both of our mod series of 
simulati ons discussed in secti on [3. 2 | as well as our previous 
work in iNelson et al. I <199St 20001. These models did not 
produce clumps, except as a consequence of insufficient res- 
olution. On t he other hand, fixed soften ing has been strongly 
advocated bv lMaver et al. I (I2002I . I20Q4T) and both their lower 
and higher resolution simulations do produce clumps. 

There are a number of questions that we must answer in 
order to understand the implications of the softening choice 
and magnitude on the outcome of a simulation. First, since 
we have no a priori knowledge of where in the disk we may 



© 0000 RAS, MNRAS 000, 000-000 



20 Andrew F. Nelson 



10 e r 




mod3TC 

Time=4246.9 yr 




mod2grv 
Time=fl82.5 yr 



n io e 



i 10 s 



10< 



1000 



100 




mod3grv 
Time=4£37.5 yr 



10 100 1 

Radius (AU) 



10 



100 



Figure 11. The cumulative maximum surface density reached at each point in the disk, binned as a function of radius in the 
mod2TC/mod3TC and mod2grv/mod3grv simulations. As in figure I1UI the azimuth averaged surface density shown in figure |S| as well 
(lower curve) and the three smooth curves show the condition from equation llll with a value of N TCBO of 20, 120 and 240 particles (top 
to bottom in each panel). 



expect fragmentation to occur, if it is to occur at all, we 
must decide how to choose a value for the gravitational soft- 
ening that allows the best reproduction of the real system. 
Given such a choice, what is the difference between the influ- 
ence that fixed softening has on the simulation compared to 
variable softening? For example, to what extent can one or 
the other choice suppress fragmentation in models where it 
should not occur, but which may be under resolved? To what 
extent are physically valid evolutionary signatures also sup- 
pressed? Can one or the other choice actually instigate frag- 
mentation in models that are otherwise stable? Can some 
values of fixed softening length both artificially suppress and 
enhance fragmentation in different parts of the same disk? 

4.2 The influence of fixed softening on 

fragmentation in simulations with insufficient 
resolution 

In this section, we investigate the influence of particular 
fixed values of the softening with a series of simulations (de- 
noted the fix TC and fix grv series in Table [3J that are iden- 
tical to the mod2TC and mod2grv simulations which we have 



shown to be unstable to numerically induced fragmentation 
with variable softening. These simulations implement a fixed 
softening using the same SPH spline kernel as is used for the 
hydrodynamic evolution, but are run, in the first case, with 
the modified kernel derivative for the pressure force calcula- 
tion or, in the second case, with the modified 2D softening. 
Each of the four pairs of simulations uses a different soft- 
ening length corresponding to the magnitude of the initial 
smoothing length at different locations in the disk. These 
are, respectively, the inner disk edge (0.5AU), the orbit ra- 
dius corresponding to the region most susceptible to clump 
formation (15AU), and near the middle of the radial span 
of the disk (30AU), and its outer edge (50AU). 

We show the mophologies of these simulations in fig- 
ures^] and EH an d the azimuth averaged surface densities 
in figures FPfl and [Tol Overall, the behavior and morphology 
observed for each variant is quite similar to that of the other. 
With small softening (fixlTC/fixlgrv) fragmentation occurs 
at essentially all orbit radii, on the same time scale as re- 
quired for the disk itself to become active. Due to this rapid 
onset, we terminated it very soon after it began, so that 
radially more distant parts of the disk simply had insuffi- 
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Figure 12. The particle distributions for the fix TC series of simulations, at the end of each simulation. 



cient time to become active. We believe that fragmentation 
would occur there too if the simulation were to be evolved 
further. Only at the inner disk edge, where the softening 
and smoothing lengths become comparable, is fragmentation 
suppressed. At all other locations, it is actually enhanced 
relative to the case of variable length softening shown in 
the mod2TC/mod2grv simulations above, where clumps only 
occurred at much larger radii. The fix2TC/fix2grv pair of 
simulations were evolved for longer in time than the corre- 
sponding mod2TC version, but ultimately they too produced 
clumps. Fragmentation occurred only in the region exterior 
to ~ 15 — 20 AU, where softening was smaller than smooth- 
ing. 

In the fix3TC/grv pair, evolved with a still larger soft- 
ening, strong spiral structures developed in the outer half 
of the disk, while comparatively weak spiral structure devel- 
oped closer to the star. Significantly, the disk was evolved 
for 5Td and did not produce any clumps over that time. For 
simulations fix^TC/fix^grv, with the largest softening, spiral 
structure is largely suppressed throughout the entire disk, 



relative to their smaller softening cousins and to mod2TC. 
The gross structure of the spiral arms that form also appear 
somewhat weaker than those in the mod2TC/mod2grv sim- 
ulations, however in the absence of a quantitative analysis 
of the pattern amplitudes (beyond the scope of this work) 
this statement remains somewhat subjective. 

From these simulations, we can make the rather unsur- 
prising conclusion that with a large enough softening value, 
clumping can be suppressed in simulations where it would 
otherwise occur due only to insufficient resolution. While 
not surprising, it is still important to quantify the both the 
value of 'large enough' and what changes in behavior occur 
as a function of gravitational softening, in order to be able 
to separate out real and artificial effects. In this case, 'large 
enough' gravitational softening appears to be larger than the 
hydrodynamic smoothing lengths in each part of the disk. 
In order to completely suppress clumping, softening must 
therefore be set comparable to the smoothing values in the 
outer part of the disk. 



© 0000 RAS, MNRAS 000, 000-000 



22 Andrew F. Nelson 

fix1grvTime=0.20Td 




Log| □ ■:■ i "i =< 1 > i 
4.5DCE-Q0 



3.6F0E+D 1 ? 
3.1LCe+[|i3 

I i.9ace+d«j 



fix2yfvTime=1.42Td 




Log; DmErty ) 



fix3grvTime=5>0Td 



fix4grv Time=5,QTd 




lose ^cisfty ) 




Figure 13. The particle distributions for the fix grv series of simulations, at the end of each simulation. 



4.3 Enhancement of fragmentation in simulations 
with sufficient resolution 

While we have seen that large enough softening can alter the 
behavior of an already numerically suspect (under resolved) 
simulation, it is also important to determine its effect on 
what might otherwise be considered a well resolved simu- 
lation. To investigate this question, we have run four sim- 
ulations, denoted mlfxTC, mlfxTC, mlfxgrv and m^fxgrv 
in table |H] with identical initial conditions and resolution 
to our mod4TC/mod4grv simulations, but with fixed grav- 
itational softening. As in the last section, the value of the 
softening is set approximately to the value of the hydrody- 
namic smoothing length at a predetermined radius in the 
disk. For the mlfxTC/grv pair of simulations, we set the 
softening equal to the smoothing at the inner edge of the 
disk (0.5 AU), so that they correspond to higher resolution 
versions of the fixl TC/grv pair. For the m^jxTC/grv pair, we 
set the softening equal to the smoothing at ~ 15 AU. These 
simulations also correspond to high resolution versions of 
the fix2TC/grv simulations. 



Figure [Tol shows the azimuth averaged and cumulative 
maximum surface densities for the two TC92 simulations. 
In contrast to the mod^TC/grv realizations, the two mlfx 
realizations produce fragments in the disk in much less than 
ITd- As for their lower resolution cousins, fixlTC/fixlgrv, 
we terminated the simulations quite soon after they began 
due to the fragmentation. We expect that had evolution pro- 
ceeded further in time, fragmentation would have occurred 
in the outermost parts of the disk as well. Consistent with 
our finding in section that higher resolution realizations 
produced fragments later than lower resolution realizations 
of the same initial condition, fragmentation occurred later 
than in the fixl variants. Fragmentation was delayed longer 
in the TC92 variant and was suppressed over a slightly larger 
portion of the inner disk. 

In both cases, fragmentation only occurred in regions 
where the fixed gravitational softening was smaller than the 
hydrodynamic smoothing, consistent with the expectation 
from the discussion in section l2~7l and in BB97 for 3D sim- 
ulations. Even though the relative size of the softening and 
smoothing lengths were the same in the low and high resolu- 
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Figure 14. Azimuth averaged and cumulative maximum surface densities (as defined for figures 151 and [TUt for the fix TC simulations, 
each defined with different, fixed gravitational softening lengths, and which are unstable to numerically induced fragmentation. The 
smooth solid curve represents critical surface density defined by the equation II II with A r re ao=6-/V nc i g i 1 . The vertical dotted line in each 
figure corresponds to the radius at which the gravitational softening and the initial hydrodynamic smoothing lengths for the particles 
are equal. 



tion realizations, clumps appeared only at somewhat larger 
orbit radii than at lower resolution. In all cases, the max- 
imum surface densities in and near the clumps exceed the 
critical value determined from equation II II and must there- 
fore be considered to be numerically induced rather than 
due to physical processes. 

Results in the m4fxTC/m4fxgrv simulations, with a 
much larger softening length, were quite different from each 
other: fragmentation did not occur with the modified ker- 
nel, but did with the modified softening, again in the outer 
part of the disk where interparticle force imbalances favored 
gravity over pressure. Although we have not shown the ear- 
lier evolutionary history of the m^fxgrv simulation, we note 
in passing that its character was markedly different than 
any other discussed in this work. High amplitude, filamen- 
tary spiral structures developed in the outer disk as early as 
~ ITd , much earlier than fragmentation occurred as defined 
by a violation of the Toomre criterion. Although the struc- 
tures developed with some regularity, they were unable to re- 
main distinct, instead becoming progressively more sheared 
out, as they disappeared into the background flow. Based 
on this behavior, we conclude that the softening length cho- 
sen for this simulation was very near the boundary between 



enhancing and suppressing fragmentation, as we saw for the 
fix2grv/fix3grv simulations at lower resolution. 

If, as we have suggested, small scale interparticle force 
imbalances enhance fragmentation when the imbalance fa- 
vors gravity, then we should expect simulations with frag- 
mentation to contain a significant number of both paired 
particles and particles at small mutual separations, but few 
or none in simulations that do not fragment. Figure 1171 
demonstrates that indeed, this supposition is true for the 
two smallest softening length variants. The number of paired 
particles rises immediately and dramatically from zero to 
over 4 x 10 4 over the short span of each simulation. Further, 
while the number of pairs appears to level out or even de- 
crease at late times, the feature is actually a consequences 
of the increasing prevalence of higher order multiple particle 
groupings, such as triples or quadruples, as demonstrated by 
the lower pair of curves in the figure, representing the pop- 
ulation of triples. 

Fieure HTn shows histograms of the full neighbor distribu- 
tions of the two larger softening variants and fieure tT9l shows 
the number of pairs as a function of time for the m4fxgrv sim- 
ulation. As we suppose, a significant population of pairs and 
near pairs (i.e. with «<0.1) develops in the modified soft- 
ening variant, which fragmented, while few do in the TC92 
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Figure 15. The same as figure H4l but for the fix grv series. 



variant. Over the course of 12Tb, only a tiny population 
of near pairs developed compared to its cousin, and we ob- 
served only ~ 5 particles to become paired, a small enough 
number to be unimportant overall. While force imbalances 
must have been present at small separations in both models, 
we conclude that the enhanced interparticle pressure forces 
in the TC92 model, coupled with the higher resolution, were 
of sufficient magnitude to eliminate the numerically induced 
fragmentation that occurred when the softening length was 
smaller. 

In an apparent contradiction of our supposition, the to- 
tal number of paired particles generated in the three simu- 
lations here, which fragmented, is smaller than in the mod4 
simulation, which did not fragment. Although we have made 
no quantitative comparison, a brief inspection indicates that 
the contradiction may be resolved by accounting for the spa- 
tial distribution of pairs and of near pairs. Because the small 
softening simulations were not evolved for a long enough 
time to become fully active, the distribution of pairs is quite 
concentrated in the regions did become active and did frag- 
ment, relative to the distribution in the mod4 realization 
where pairs are found spread much more evenly throughout 
the system. The m^xgrv simulation, which did become fully 
active, generated many fewer pairs but did generate a large 
population of near pairs, indicating again the consequences 
of force imbalances favoring gravity. Given that the effec- 
tive resolution in the mod4 run is decreasing with time, we 



expect that if it were to be evolved further, and the concen- 
tration of paired particles increased further, fragmentation 
would occur. 

Disks realized with small fix softening fragmented in 
both the TC92 and the modified softening variants, and both 
in their low resolution realizations and at resolution high 
enough to satisfy the Toomre condition defined in section 
13.31 An immediate conclusion from the similarities might 
be that, since the fragmentation occurred each realization, 
the simulations are in fact converged. While this conclusion 
may in fact be correct, it would also be seriously misleading. 
Due to flaws in their design, the simulations converge to 
an incorrect result. Interpretation of that result in terms 
of the physical behavior of the system and the importance 
of other correctly implemented physical processes becomes 
difficult or impossible. While the result and its invalidity 
may be clear in cases like our test simulations, the same 
statement may not be true in other simulations where the 
results are more difficult to verify. As a trivial example, the 
high and low resolution variants with larger fixed softening 
for which both low resolution variants fragmented, but only 
one at higher resolution. Our conclusion in this case must be 
limited only to the fact that the TC92 kernel modification 
is more likely to suppress fragmentation than the softening 
modification. 

Two other, more technical aspects of simulations with 
fixed gravitational softening deserve mention. First, the cpu 
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Figure 16. Azimuth averaged and cumulative maximum surface densities for the two simulations using the TC kernel derivative (left) 
and the two with 2D softening variant (right). In each case the gravitational softening length is fixed to the smoothing length at the inner 
edge (top) of the disk or at 15 AU (bottom), and all realizations are expected to be stable according to the Toomre criterion defined by 
equation 1111 with N TeBO =6N ne [ g ^. The vertical dotted lines correspond to the radius at which the gravitational softening and the initial 
smoothing lengths for the particles are equal. 



time required to run simulations with large, fixed softening 
to some specified time is far longer than with either the vari- 
able softening case or the cases with small, fixed softening. 
The reason is that with a large softening length, the inter- 
particle spacing requires that a much larger fraction of the 
total number of gravitational interactions between particles 
be calculated as 'atoms' rather than grouped together as 
nodes, as is ordinarily done in tree-based gravity solvers in 
general use for particle simulations. Secondly, in cases where 
fragmentation does occur, a simulation with large fixed soft- 
ening can proceed for much longer in time, because the time 
step size does not decrease nearly as much as occurs in a 
variable softening simulation at the same resolution. Con- 
tinuing a simulation for a long time after clumping occurred 
and the resolution criterion has been violated would be of 
limited utility however, since the clump may cause large per- 
turbations to the rest of the system that would not occur 
otherwise. 

4.4 Recommendations for choice of softening 

We conclude from our models that a fixed softening may ei- 
ther enhance or suppress density inhomogeneities depending 
on the specific choice of softening value relative to the hydro- 



dynamic smoothing. No single value of softening may be con- 
sidered 'optimal' in a disk with a large radial extent. In the 
extreme case, an incorrect (too small) softening value may 
induce fragmentation to occur in a simulation that would 
not, given another fixed choice or a variable softening set 
to the hydrodynamic smoothing length. Given an incorrect 
but too large value, suppression of physically correct struc- 
ture may occur. In contrast, while a variable gravitational 
softening length may produce some violation of energy con- 
servation, in the examples shown here the outcomes of the 
simulations are not affected in a manner as drastically as 
with fixed softening. We therefore recommend that simu- 
lations involving self gravitating hydrodynamic systems in- 
corporate a spatially and temporally variable gravitational 
softening, whose length scale is the same as that over which 
the hydrodynamic quantities are smoothed. 

Our recommendation is similar to, but stronger than 
that made by BB97 because we observe that fragmentation 
may develop in simulations whether or not the Toomre cri- 
terion (eauationlolor tTni is satisfied. Rather than signaling a 
temporal endpoint, beyond which BB97 conclude that evolu- 
tion towards some final collapsed state is under resolved, we 
find that force imbalances at small particle separations can 
alter the evolutionary trajectory of a system which would 
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Figure 10. The cumulative maximum surface density reached 
at each point in the disk, binned as a function of radius in the 
modS and mod4 simulations. For reference, the azimuth averaged 
surface density shown in figure|3as well (lower curve). The three 
smooth curves show the condition from equation II II with a value 
of iVreso of 20, 120 and 240 particles (top to bottom in each panel). 
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Figure 18. Histograms of the neighbor separations for all neigh- 
bors of all particles in the m4fxTC (top) and m4fxgrv (bottom) 
simulations. 
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Figure 17. The number of multiple particles in the mlfxgrv 
(solid) and mlfxTC (dotted) simulations as a function of time. 
The two curves ending near 4 X 10 4 represent the number of parti- 
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Figure 19. The number of particle pairs present in the m4fxgrv 
simulation as a function of time. 



otherwise never approach gravitational instability into one 
that does. Given that BB97 did not specifically study sys- 
tems starting in a dynamically stable equilibrium condition 
for which fragmentation is known not to occur, as found in 
our disks, the behavior we observed in our simulations may 
not have been observable in their work. We suggest that it 
is independent of the specific system being simulated, and 
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in fact applies to the Jeans collapse situat ion as w e ll. Ou r 
conclusion is also consistent with that of iDehnenl i200lt) . 
who made a detailed study of the effects of various choices 
of softening on the force errors produced by a static distri- 
bution of particles, and who recommends spatially variable 
softening in that case as well. Since in nearly all particle 
based hydrodynamic simulations, however, the smoothing 
length is not fixed in time, our recommendation goes be- 
yond Dehnen's study of a static distribution to including 
temporal variation as well. 

We base our recommendation on the results of our sim- 
ulations and on the following argument. Completely inde- 
pendent of whether or not either the gravitational or pres- 
sure force actually can be resolved on better or worse length 
scales than the other, is the consequence of the assumption 
underlying both softening and smoothing. Namely, an error 
in the net force is present by assumption when the contri- 
bution from either source is not fully resolved. We there- 
fore have only the choice of how to handle this error most 
gracefully in the code. We submit that avoidance of large 
force imbalances in regions that are by definition under re- 
solved is a trait much to be desired in the properties of the 
code because pathological behavior that results from such 
imbalances-artificial fragmentation or suppression of physi- 
cally realistic fragmentation-are avoided. 

Such a recommendation is not without cost however 
because it is equivalent to recommending that simulations 
using SPH should not require that energy be conserved. Al- 
though we have not attempted to quantify the magnitude of 
the violation, in general we believe it to be a relatively small 
component of the total energy of a given particle since only 
neighbor particles will contribute. Also, the violation will 
be smaller in higher resolution simulations because those 
neighbors will represent a proportionally smaller component 
of the total mass in the system, from which the potential 
energy of a given particle is derived. Finally, the fact that 
the gravitational forces and potential are approximated by a 
summation of terms derived from a tree search already vio- 
lates strict energy conservation. Given the large body of lit- 
erature using such approximations, often checked for verac- 
ity against other methods or observations, it would appear 
that few, if any, unacceptable consequences originate from a 
low level of violation. Assertions such as these are no replace- 
ment for a m easurement however, and we refer the reader to 
iBenzl Jl990l) for discussion of a method to make an approxi- 
mate quantification. We are also aware of current work, de- 
velopi ng on the techniques discussed in lPrice fc Monaehanl 
(2004) and by the same authors (Price and Monaghan 2006, 
submitted), to create an SPH formulation that allows vari- 
able softening, while still conserving energy. Adapting SPH 
codes in existence to utilize this technique would appear to 
be highly desirable, since it would remove the most serious 
drawback of variable softening. 

We have considered and discarded a third, intermedi- 
ate softening option, in which the gravitational softening of 
each particle is set, for example, to the initial value of the 
hydrodynamic smoothing length, and is thereafter fixed. Al- 
though no specific tests have been performed to study this 
question, we believe that using such a prescription will yield 
results similar to the fixed softening case with a small soft- 
ening length. For example, a particle initially located in a 
dense region (small softening value) , which later moves into 



a less dense region (larger softening value) could act as an 
attractor for other particles in its current neighborhood. An 
exactly opposite effect might be seen for a particle initially 
in a low density region which moves to a high density region. 

Secondarily, we question the aesthetic utility of such 
a prescription because it implies knowledge of the internal 
mass distribution of each particle that is both independent 
of any hydrodynamic expansion or contraction of the gas 
(the assumption underlying a temporally variable softening) 
and needlessly complex. The complexity is due to the fact 
that using constant, individual softening lengths suggests 
that information is known about not only the internal mass 
distribution of particles taken as an ensemble, but also the 
information that each particle has a different internal mass 
distribution from every other. We can postulate no circum- 
stances where such exquisitely detailed information is likely 
to be available. 



5 SIMULATIONS IN 3 DIMENSIONS AND/OR 
USING GRID METHODS 

The wave analysis underlying our resolution criterion makes 
the assumption that the flow is two dimensional: an inte- 
gration over the z coordinate is assumed. The test problem 
defined in section El w as chosen to be simil ar to the models 
originally employed bv lNelson" et al. Hl99gfi . and was there- 
fore performed in the same two dimensional approximation 
as the originals. This approximation has historically been 
common in other works as well. In this section, we examine 
methods for applying the criterion to fully three dimensional 
models, and to the question of whether 3D simulations re- 
quire additional constraints on resolution in order to ensure 
their accuracy. 

5.1 Calculating surface density in 3D simulations 

For fully three dimensional simulations, where a surface den- 
sity is not directly available, the Toomre criterion cannot 
be applied without some additional computation. For grid 
based simulations, the computations are trivial and involve 
only an integral of the volume density in the z direction. In 
cylindrical or Cartesian coordinates, for example the inte- 
gral would pass to a summation of pkSz over a set of grid 
zones with spacing Sz, with pk being the volume density in 
the fcth grid zone in the z direction. 

For 3D SPH simulations, the computation of a surface 
density is more complex because there is no grid. Instead, we 
recommend the creation of a temporary grid onto which the 
volume density can be mapped. Then the surface density 
can be computed as a summation in the z direction as in 
the case for grid based simulations. In test models, we have 
found good results by using a mapping to the grid based on 
the same smoothing technique used to determine the volume 
densities at the particle locations, for the mapping. At the 
location of each grid zone, we determine a list of 'neighbor' 
particles which may contribute to the density of that zone, 
then we use the formula 

p(i,j,k)= m P W(r,h) (26) 

neighbors 

where m p is the mass of each neighbor particle, W is the 
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smoothing kernel, r = |r p — r zonc | is the distance between 
the particle and the grid zone at location k) and h is the 
mutual smoothing length of the particle and zone, defined 
as max(ft ZO nc, h p ). The grid smoothing length is defined as 

/izonc = y/5r? + (r i 8<l> j ) 2 +8z% (27) 

where the components are the grid spacing in each of r, 
<j> and z. We include a smoothing length for the grid itself 
in order to ensure that every particle contributes to at least 
one grid zone even in the case of vanishingly small smoothing 
lengths. 

For the test problem discussed in sections 15.31 15.41 and 
!5.6l below. we map the volume density onto a cylindrical co- 
ordinates grid where the radial coordinate is logarithmically 
spaced and the azimuth coordinate is uniformly spaced so 
that its spacing is 8r « r8<f> everywhere. In order to resolve 
the smallest scale features of the flow, we have chosen our 
grid resolution so that its spacing is comparable to the small- 
est smoothing lengths in the calculation. This grid spacing 
will in practice over resolve the features of most of the flow, 
and will not strictly be required in most circumstances. Our 
choice is motivated by the desire to ensure that the choice 
of grid (i.e. h zonc ) will contribute negligibly to the derived 
densities and the quality of the fitted quantities in our anal- 
ysis that depend on them. In any case, the computational 
expense of choosing of a finely spaced grid for the mapping 
is small in terms of the time to run a simulation because the 
mapping is not required for the simulation itself, but only 
for post-processing. 

A conceptually simpler mapping would be to use a Par- 
ticle in Cell technique to assign the mass of each particle to 
the grid cell that contains it, thus determining a volume den- 
sity for that cell. Unfortunately, we cannot recommend this 
method because it is susceptible to large errors due to the 
fact that SPH particles may have smoothing lengths much 
different than the grid spacing. For example, when a par- 
ticle's smoothing length is larger than the grid spacing, it 
should contribute to the density over a larger region than 
assumed by the PIC technique. 



5.2 An unwarranted approximation: p — T./2H 

When two dimensional simulations are performed, it is com- 
mon to obtain an approximation to the volume density by 
dividing from the disk's surface density by some factor mul- 
tiplying the disk's scale height. For an non-self gravitating 
isothermal gas, the factor will be V2ty « 2.5 (see section IQ1 
below), but will vary slightly depending on the exact input 
physics. For purposes of this section and for simplicity, we 
shall assume a factor of two exactly. As discussed at the end 
of section 12.41 the Jeans or Toomre stability criteria may 
be applied using this conversion to obtain one or the other 
of surface or volume density as required. In this section, we 
will investigate the impact that this conversion may have on 
the indirectly obtained values used in the Jeans criterion. 
We will refer here to the criteria directly applied through 
the use of equations [5] or [8] as the Jeans or Toomre criteria, 
while the criteria applied through the use of the disk scale 
height to convert surface to volume density (or vice versa) 
as the approximate Jeans or Toomre criteria. 

In order to proceed, it is necessary to examine the re- 



sults of a simulation modeled fully in 3D, so that both vol- 
ume density and surface density are known. In one of a s eries 
of papers studying fragme ntation in d isks, iBossI i2002T) dis- 
cusses 3D disk models and IBossI i2004j) continues this study 
with a close examination of the properties of a few models 
(primarily the model named 'edh') from the earlier work. 
We will use these models as a test bed for our investigation. 
Since our study has not specifically defined a value for T in 
equation |S] we will assume that its value is T = 1/4, which 
is identical to the value of J required for Jeans collapse sim- 
ulations. While we believe this to be a reliable assumption, 
only by examining the results of a specifically chosen test 
problem can we be fully confident that the value is suffi- 
cient. The consequences of the uncertainty in specifying T 
are that comparisons between the magnitudes of T and J 
and the critical lengths associated with them, may be more 
difficult to interpret. In this section however, we seek pri- 
marily to test the veracity of the approximate forms of the 
criteria compared to the exact forms, and the absolute mag- 
nitudes are less important. 

Boss's model consists of a circumstellar disk modeled 
between 4 and 20 AU, evolved on a grid of 100 x 23 x 256 (or 
x512 in some simulations) zones asymmetrically distributed 
primarily near the disk midplane. He includes an ideal gas 
equation of state and radiative cooling in the diffusion ap- 
proximation, but omits an artificial viscosity, thereby omit- 
ting the viscous and shock heating modeled by it. With this 
model, he evolves a disk for 345 yr, at which time the sim- 
ulation begins to violate the Jeans criterion due to onset of 
clump formation, and is terminated. We will apply our sta- 
bility criteria to this model using each of the four relevant 
versions discussed in section 1231 

The top panel of figure |2T]1 shows the value of the criti- 
cal wavelengths at the time the simulation was terminated, 
using each of the four applicable implementations of the sta- 
bility criteria. The failure of the Jeans criterion was used 
to define the end of the simulation, so its minimum value 
matches the grid spacing very closely. The exact Toomre cri- 
terion yields a critical wavelength that is about 30% smaller 
than that from the Jeans analysis. Using it and assuming 
that T = J, the model has already passed the boundary 
beyond which it becomes numerically suspect. Due to the 
speed at which clumping occurs, once underway, we do not 
believe that this small difference is very significant however, 
except to note that the simulation would have been termi- 
nated slightly earlier and at a lower density. At all other 
locations, the Toomre wavelength exceeds the Jeans wave- 
length and the evolution of the Boss simulation leading up 
to its termination passes the resolution test. The conclusions 
made from it stand on valid numerical grounds according to 
the Toomre resolution criterion. 

On the other hand, the approximations of the Jeans 
(light solid) and Toomre (light dashed) criteria lead to criti- 
cal wavelengths in the forming clump that differ by a factor 
of 1.56 larger than and 4.6 smaller than the grid spacing, 
for the approximate Jeans and approximate Toomre crite- 
ria, respectively. Moreover, using the approximate Toomre 
criterion, the simulation is already well past the time for 
which its evolution is numerically valid. Each of the approx- 
imate forms uses the isothermal scale height to render the 
conversion. Does this quantity remain relevant in a collaps- 
ing region? 
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Figure 20. (a) Critical wavelengths for the Boss (2002) sim- 
ulation along the radial column including the forming clump. 
The exact Jeans and Toomre criteria are shown with heavy solid 
and dashed lines respectively, while the approximate Jeans and 
Toomre criteria (as defined in the text) are shown with light solid 
and dashed curves, respectively. The horizontal solid line, defines 
the spacing of the grid on which the system was evolved, (b) The 
ratio of the Jeans and Toomre wavelengths to the isothermal scale 
height (H = c s /fi) as functions of orbit radius. 



The bottom panel of figure |2Tj1 shows the ratios of the 
Jeans and Toomre wavelengths to the scale height. Except 
near the forming clump, the Jeans wavelength is many times 
the disk scale height, making its application as a resolution 
criterion unclear according to the argument in section 12.31 
The Toomre wavelength is also many times the disk scale 
height, but this only makes its application to a resolution cri- 
terion more appropriate, according to the same arguments. 
We therefore recommend that it be used in favor of the Jeans 
criterion in simulations of disks. In and near the collapsing 
region, both wavelengths become significantly smaller than 
the scale height. Since those wavelengths loosely define the 
physical size of the collapsing body, we can conclude that 
the structure of that object is no longer well represented as 
a slightly perturbed disk structure under which assumption 
the scale height was originally defined. Thus, the approxi- 
mate criteria should not be used because their values be- 
come inaccurate near the clump, due to the fact that the 
disk vertical profile becomes distorted. 



A simple 3D test problem to test the accuracy 
of the hydrodynamics in disk simulations 



5.3 



In section 12781 we discussed the possibility that a disk simu- 
lation performed in 3D may not in fact be truly three dimen- 
sional because the disk thickness may be comparable to the 
smoothing lengths of the particles in and SPH simulation, 
or the grid dimension in a grid based simulation. Here we 
explore the consequences of such simulations using a simple 
test problem. Specifically, the problem of a non-self gravi- 
tating disk evolved with an isothermal equation of state. For 
this system, an analytic expression for the vertical structure 
of the disk can be derived exactly, and direct a comparison 
between the analytic result and theory can be made. This 
test problem is equally applicable to grid based simulations 
and to particle simulations, and we propose it as a general 
test for interested numericists. 

The vertical structure of the gas in a non-self gravitating 
disk will obey the equation 



dp _ GM t z 

Tz ~ P ( r 2 + 2 2)3/2 



(28) 



where p and p are the pressure and volume density, z is the 
altitude above the disk midplane and r is the cylindrical 
radius. For an isothermal equation of state, the pressure is 
related to the density and sound speed through the definition 
p — p(? s . Using the equation of state to replace pressure with 
density, the solution to the differential equation is 



p(z) = p e" 



(29) 



where H = c a /fl defines the isothermal scale height and 
n = \J GMt/r 3 . For any given location in the disk, the 
surface density can be obtained by integrating equation 1291 
over all z, giving the relation 



Po = 



2-kH 



(30) 



We use a modification of the 2D test problem described 
in section mi for this test. Particle layout, temperature and 
surface density are defined as in the 2D problem, such that 
a minimum Toomre Q value of ~ 1.5 would be obtained if 
self gravity were included. To obtain an initial z coordinate 
for each particle, we use a Gaussian pseudo random number 
generator whose width is defined by the local value of the 
disk scale height, consistent with equation 1291 Definition of 
the initial velocities takes account only of stellar gravity and 
pressure forces rather than disk self gravity as well. Veloci- 
ties in the z coordinate are set to zero. 

Before proceeding further, it is significant to note that 
the specific model we propose is a step removed from those 
of greatest relevance for actual systems because vertical den- 
sity distributions depend on the details of the physical model 
employed. Also, as a practical matter, the exact definition 
of 'scale height' itself can lose meaning in such models. To 
the extent that it does remain meaningful, any resolution re- 
quirement may be sensitive to the details of the model if, for 
example, a larger fraction of the mass were in an extended 
envelope or were concentrated closer to the midplane than 
with an isothermal structure, as is the case for the isentropic 
structure discussed in section 15.71 below. In the context of 
the spatial distribution of particles or grid cells that repre- 
sent the mass and at a coarse level of examination, models 
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that result in a higher midplane mass concentration will be 
somewhat similar in structure to colder isothermal disks in 
which the scale height is simply not as large. We therefore 
believe that any modifications to the critical resolution ap- 
propriate for other models will be not be large and proceed 
to define a criterion for isothermal structure that we believe 
will be generally applicable. In the following two sections, we 
will however, proceed to define a resolution criterion in two 
complimentary ways: as a requirement per scale height at 
the disk midplane, and as a resolution per vertical column, 
in order to facilitate its use more generally. 



5.4 Results from simulating our 3D test problem 
with SPH 

We have run a series of simulations modeling the 3D test 
problem at different resolutions, from 1.3 x 10 5 to 10 6 parti- 
cles, each labeled sgoff3-6 and defined in tabled We allow 
each simulation to run for two orbits of the outer edge of the 
disk in order to ensure that the configuration represents the 
true system as realized by the hydrodynamic code rather 
than any peculiarity of the initial configuration. We then 
map the volume density and velocities onto a three dimen- 
sional grid and calculate the surface density at each location. 
We pose the following question for the configuration, which 
must be answered affirmatively before the results can be val- 
idated: does the configuration generated by the simulation 
correctly reproduce the correct midplane densities and scale 
heights everywhere in the disk? 

To answer the question, we generate least squares fits 
to the densities as mapped onto a temporary grid to com- 
pare to the analytical values defined by equation 1291 Each 
fit returns three parameters corresponding to the values of 
po and H in equation 1291 We limit the fits to the radial 
region between 1 AU and 45 AU in order to eliminate the 
possibility that conditions at the inner or outer boundary of 
the disk distort the true picture of the structure. Inspection 
of the disk structure shows that these limits ensure that the 
fits are limited to the regions where the rotation curve is un- 
affected by either the gravitational softening of the star at 
the inner boundary or large density gradients near the outer 
boundary. In order to remove any potential systematic er- 
rors derived from the resolution of the grid or simulation 
itself, we average the volume densities at each altitude over 
patches of approximately one scale height on a side, so that 
the same number of fits are generated for each simulation. 
Finally, in order to retain self consistency in the face of any 
evolution of the disk away from the initial surface density 
profile, we derive the analytic values for po and H using the 
locally determined value of the surface density and rotation 
velocity, through the relation H = c 3 /Q, and equation 1301 

Patches are defined in successive radial rings of the disk 
by the condition that a ring's radial extent is the local, an- 
alytically determined value of the scale height. Their az- 
imuthal extent is determined by the same condition, so that 
the total number of patches in a ring is defined by its circum- 
ference divided by the scale height. At the grid resolution 
chosen for the analysis, each patch covers a region of ~ 200 
or more cells in the vertical coordinate, and 10 x 10 radial 
and azimuthal grid cells near the inner disk edge where the 
scale heights are small, and 15 x 20 near its outer edge where 
they are much larger. In total, approximately 14000 patches 
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Figure 21. Histograms of the best fit midplane density, po (top) 
and scale height, H (bottom). The black, red, green and blue 
curves correspond respectively to the progressively higher resolu- 
tion simulations sgoff3-sgoff6. The bins of the histogram are of 
width 0.25%. 



are required to cover entire disk, so that 14000 separate fits 
of the vertical structure in the disk are performed. 

Figure 1211 shows histograms of the ratios of the fitted 
quantities po and H to their analytic predictions, normal- 
ized to the total number of patches. At the lowest resolu- 
tion shown, the vertical structure in the simulation is clearly 
more extended vertically than predicted, with the distribu- 
tion of scale heights extending more than 50% larger than 
the analytic values and the midplane density distribution 
extending nearly 50% below the analytic values. Both dis- 
tributions are spread relatively evenly over the range ex- 
tending from the most extreme under or over estimates on 
one end, to nearly the correct analytic value on the other. 
The distributions become narrower for each of the each of 
the higher resolution variants of the model, and appear in 
fact to be converging to the analytically expected values. 
However, even at a resolution of one million particles, the 
highest resolution in our study, the peaks of the distribu- 
tions deviate from the analytical values by about 5% and a 
significant population of patches are best fit with midplane 
densities and scale heights as much as 20% below or above 
the analytic values, respectively. 

We can convert the correspondence between simulation 
and theory into a resolution requirement if we can quantify 
the minimum number of particles required in a vertical col- 
umn for which the fit parameters are accurately reproduced. 
The top panels of figure l^ show the vertical resolution of the 
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Fit/Analytic Midplane Density Ratio 

Figure 22. Top panels: Scale height to smoothing length ratio for all particles in the sgoff4 and sgoff6 simulations, as labeled. Bottom 
panels: Histograms of the ratio of the best fit to analytic midplane densities expected for all patches in the disks between 1 and 45 AU. 
The color of each histogram corresponds to to patches in the ranges 1-3 AU (black), 3—9 AU (red), 9—27 AU (green) and 27-45 AU 
(dark blue), as shown in the top panels. Patches at orbit radii smaller than 1 AU or larger than 45 AU (light blue) are excluded from 
the histograms in order to minimize distortions due to decreased resolution at the inner and outer disk edges. As in figure |2T1 histogram 
bin widths are set to 0.25%. 



disks quantified as a ratio between the expected scale height 
and the smoothing length of each particle. The bottom pan- 
els show the quality of the reproduction of the analytically 
determined disk midplane densities as in figure |2"TI broken 
out into separate histograms for each of four radial regions 
of the disk. 

The vertical resolution varies between less than one 
smoothing length per scale height in the inner disk to as 
many as six near the outer edge of the high resolution sim- 
ulation, sgoff6. The ratio takes values over a range from 
less than unity up to a orbit radius dependent maximum 
value because particles located near the midplane at some 
radius, where densities are high, will have correspondingly 
smaller smoothing lengths. Particles located at higher alti- 
tudes where densities are lower will have larger smoothing 
lengths. Dividing each by the scale height (constant at a 
given radius) produces the range. The maximum value of 
the ratio varies as a function of orbit radius primarily be- 
cause the scale height itself is a function of radius, through 



the combination of the predefined temperature profile (equa- 
tion and the rotation curve. 

In the 1-3 AU region of the disk, the entire vertical 
structure of the disk is represented by the smoothing length 
of only a single particle. The exact positions of neighboring 
particles relative to each other therefore cause a correspond- 
ingly larger influence on the densities, as demonstrated by 
the very wide distribution of midplane densities present in 
both the sgoff4 and the sgoff6 realizations, relative to their 
analytic values. In each successively more distant region, and 
at both resolutions, the distribution becomes narrower and 
of better quality, demonstrating the value of the increased 
vertical resolution. For example, the vertical structure is re- 
solved by ~ 1 — 2 smoothing lengths per scale height at 
r < 3 AU, and the fitted midplane densities fall between 65 
and 95 percent of their analytic value for the 260000 particle 
simulation and 75 and 100 percent in the 1 million particle 
run. At the other end of the spectrum, density fits fall within 
~ 5% of their predicted values only in the two outermost re- 
gions in the high resolution simulations, corresponding to 
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orbit radii > 9 AU (i.e. the green and dark blue histograms 
in the figure) , and only part of the distribution in the single 
outermost region of its lower resolution counterpart, corre- 
sponding to radii >27 AU. 

If we specify that midplane densities within 5% of their 
analytic values are sufficiently accurate for the purpose of 
simulating the evolution of a circumstellar disk, then only 
these outer regions of the disk have sufficient accuracy. In 
the two best resolved regions of simulation sgoff4, ratios of 
fit to analytic midplane density values peak near ~ 94 — 95% 
and extend to as low as 90%, resolving one scale height with 
at most ~ 3 — 4 particles. In the same regions of sgoff6, The 
peaks of the distributions increase to near ~ 96 — 97%, with 
tails extending to as low as ~ 92 — 93%, resolving one scale 
height with as few as 3.5 particles at 9 AU and as many as 
6 further out. We conclude that midplane densities within 
~ 5% of their analytic values are achieved only in regions 
where the vertical structure of the disk is resolved with at 
least ~ 4 smoothing lengths per scale height at the disk 
midplane. 



5.5 Interpreting the vertical resolution 

requirement at the disk midplane in the 
context of the full thickness of the disk 

It is important to note that the requirement for > 4 smooth- 
ing lengths per scale height near the disk midplane is suffi- 
cient only for obtaining reasonably accurate midplane densi- 
ties, and only for a mass distribution described by equation 
1291 Quantities at higher altitudes may still be significantly 
under resolved because smoothing lengths are larger there 
due to the correspondingly lower densities. Errors due to 
insufficient resolution at high altitudes will be of lesser im- 
portance for simulations modeling gravitational instabilities 
because the magnitudes of the pressure and gravitational 
forces important for fragmentation will be largest where the 
mass is concentrated, close to the midplane. Nevertheless, 
some particles must be present at high altitudes in order to 
compress those lower down to the required densities. How 
many are enough? 

In this section we quantify the minimum number of par- 
ticles required per vertical column needed to obtain good 
correspondence with theoretical expectations for density 
near the midplane in our models. Because it is a metric 
that is effectively an integral over a full vertical column, 
this quantification is likely to be much less sensitive to the 
exact distribution of mass over the column, and so more 
generally applicable than a quantification per scale height. 
It also eliminates the possibility of misperceptions based on 
the idea that the total mass outside the midplane is small 
enough to be neglected in estimates of the total resolution 
required per vertical column. 

For the isothermal disks in this study, 98% of the disk 
mass will reside within three scale heights of the midplane, 
according to equation 1291 meaning that the mass is effec- 
tively distributed over a total of six scale heights accounting 
for symmetry above and below the midplane. At the simplest 
level, we might therefore approximate the required vertical 
resolution to be ~ 24 smoothing lengths per vertical column 
to avoid serious errors in the midplane densities. A some- 
what more accurate quantification might also account for 



the fact that mass is preferentially located near the mid- 
plane, with fewer particles present at high altitudes. 

In principle, we can derive an approximate number of 
particles required per vertical column if we simply quantify 
the particle density over all altitudes in a single vertical col- 
umn. At the minimal resolution required to obtain accurate 
midplane densities however, we must presume that the high 
altitudes remain inaccurate due to the fact that densities are 
much lower there and particle separations larger, and they 
may therefore remain under resolved. Any direct quantifica- 
tion of the resolution required based on the actual distribu- 
tion of particles in a simulation would therefore suffer from 
the inaccuracies in their high altitude distribution. 

A fact that allows the analysis to proceed is that, for the 
purpose of ensuring accurate midplane densities, only the 
total weight of particles at high altitudes is needed rather 
than their distribution. This is important because weight 
is a quantity integrated over the vertical column of mate- 
rial and will therefore take the same value whether or not 
the mass distribution responsible for it is well resolved. Any 
high altitude mass distribution yielding the same weight will 
yield the same midplane density so that, given a correct mid- 
plane density, we may infer that the total weight is approx- 
imately correct as well. We may therefore assume that the 
high altitude mass distribution derived from our simulations 
is correct (even though it may not be), for the purpose of 
converting densities at a given altitude to a corresponding 
number of particles needed to supply the correct compres- 
sional forces to midplane material. 

We can quantify the total particle count per vertical 
column in the simulations using a modification of the patch 
averaging strategy used above. Rather than creating full 3D 
density distributions for each patch, it is sufficient only to 
quantify the averaged number of particles located in each, 
and to create a coarse vertical distribution of the particles 
by assigning each particle to a histogram bin according to its 
altitude. The result is the patch averaged volume density of 
particles in each histogram bin, from which a linear (vertical) 
particle density per scale height can easily be obtained by 
taking a cube root. A final link is to note that interparticle 
separations are ~ h and that, by definition, the histogram 
bins are of width H, so that the linear particle density is 
equivalent to the ratio of smoothing length to scale height. 

Figure l551 shows histograms of the patch averaged ratios 
of H/h defined using this procedure, separated into the same 
radial regions defined in section |5~4*1 Consistent with expec- 
tations based on the physical model and on the fit results 
above, each histogram displays a maximum at the midplane 
and non-zero populations to as high as ~ 3 — 4 scale heights 
in each direction. Histograms with larger net populations 
exhibit a progressively more peaked structure, reflecting the 
better agreement with the analytical model provided by the 
higher resolution. For the present purposes, it is interesting 
to note that in every region of the disk, the total resolution 
per vertical column is several times the resolution assuming 
that only the midplane contributes significantly to the total. 

In section 15^4*1 we determined that only the two outer- 
most regions of the sgoff6 simulation were well resolved in 
our simulations. Here, we see that on average, the resolution 
of vertical columns in these regions is ~ 17 and ~ 22 par- 
ticles per column for the 9-27 AU and 27-45 AU regions, 
respectively. For the lower resolution sgoff4 simulation, the 
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Figure 23. The patch averaged vertical, linear particle den- 
sity binned in histograms of width H, for the sgoff4 and sgooff6 
simulations. From top to bottom, the histograms in each panel 
correspond to the dark blue (solid), green (dotted), red (short 
dashed) and black (long dashed) regions as defined in figure |2*21 
The numerical values associated with each curve are sums over 
all histogram bins thus defining the averaged total particle counts 
per vertical column. 



highest overall vertical resolution reaches ~ 14 particles per 
column, and all other regions in both simulations fall pro- 
gressively farther below this value. Since fits in all of these 
regions deviate from the analytic expectation by more than 
5%, we therefore conclude that a minimum of ~ 17 — 20 
particles per vertical column are required to adequately re- 
produce the density structure at the disk midplane. 



5.6 Significance of the results of the 3D test 
problem 

In self gravitating systems, the exact balance of pressure 
and gravitational forces will tip that system towards or away 
from two dramatically different outcomes: fragmentation or 
continued smooth evolution. In marginally stable systems 
of interest to modelers, the exact balance between the two 
forces will be determined from terms of nearly equal mag- 
nitude, but of opposite sign. If one of those quantities is 
incorrectly calculated, the outcome will be dramatically dif- 
ferent. 

The results of the test above demonstrate that simula- 
tions with too few particles tend to underestimate the mid- 



plane densities, in the case of 2.6 x 10 5 particles, by as much 
as 30-35%. This is important because the pressure will be 
underestimated by a similar factor through the equation of 
state. Fragmentation will therefore be enhanced in an un- 
der resolved simulation compared to that of either a well 
resolved simulation or, more importantly, a real, physical 
system. Because we have not performed self gravitating sim- 
ulations, we cannot specify the exact effects the errors will 
have in simulations or the level of enhanced fragmentation 
that may occur however. Such details are difficult to specify 
with precision because any enhancements may be mitigated 
in part by the fact that the cause of erroneously low densi- 
ties in SPH is simply that particles separations are greater 
than they should be, so that gravitational attraction be- 
tween them is correspondingly less. Any possibility of miti- 
gation of this sort provides little comfort relative to simply 
solving the hydrodynamic equations accurately in the first 
place however. 

For the disk morphology discussed here, resolution suf- 
ficient to accurately reproduce the disk's vertical structure 
and midplane density beyond ~ 10 AU requires at least one 
million particles. Although we have not attempted to quan- 
tify the differences through simulations, colder disks closer 
to the fragmentation boundary will have correspondingly 
smaller scale heights, and will therefore require still more 
total particles to ensure adequate resolution. 

For example, the disk scale height and Q values both 
depend directly on the sound speed, so a change in one will 
be reflected proportionally in the other. A disk with min- 
imum Q = 1.3orQ = l.l will decrease the scale height 
below that in our simulations by a factor 1.3/1.5 ~ 0.87 
or 1.1/1.5 ~ 0.73 respectively. In order to retain adequate 
resolution of the vertical structure, smoothing lengths must 
be decreased by a similar factor by increasing resolution. 
In 3D, the magnitude of the increase will be a factors of 
1/0. 87 3 » 1.5 or 1/0.73 3 w 2.5 corresponding to ~ 1.5 or 
~ 2.5 million particles. No simulations of circumstellar disks 
in the context of planet formation have yet been performed 
at such high resolution. 

Even at a resolution of 1 million particles, the vertical 
structure in our simulations is not accurately reproduced in- 
side 10 AU, where midplane densities drop well below their 
analytic values. This is important for simulations of gravita- 
tional fragmentation because the Jovian planets in our own 
solar system formed at such radii. The erroneously low den- 
sities may lead to enhanced fragmentation in those regions, 
even if the structure is accurately modeled further out. Also, 
the resolved and unresolved parts of the disk are in no way 
isolated from each other. As waves or other structures prop- 
agate through the disk, their subsequent evolution will be 
perturbed by the change in resolution, so that the evolution 
throughout the entire simulation become suspect. 

Both the small values of the scale height at small orbit 
radii and the insufficient resolution there are consequences of 
the r -1 / 2 temperature profile assumed for the disk, and com- 
monly used throughout the literature. Through the sound 
speed and rotation curve, the scale height will exhibit a pro- 
portionality H/r oc r 1//4 and, for minimum Q — 1.5, will 
take values of ~ 0.025 near 1 AU increasing to ~ 0.065 
near 45 AU. The temperature profiles in our simulations 
reflect the stage of active disk evolution during which grav- 
itational instabilities are most likely to be present. A tem- 
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perature profile inversely proportional to radius will yield a 
flat profile and so avoid the issues of variable vertical resolu- 
tion, but such steep pr ofiles are not observed in real systems 
jBeckwith et al. Ill99d) . and so we discount them here. 



5.7 Additional requirements on vertical resolution 
for simulations including radiative transfer 

Above, we showed that the vertical structure of an accre- 
tion disk must be resolved with some minimum number of 
particles if mass densities in the disk midplane are not to 
be substantially underestimated. Although this requirement 
may be sufficient to ensure accurate simulations in the case 
when comparatively simple physical models are employed, 
it will be less so when radiative transfer is included. Here 
we describe an exercise meant to illustrate requirements for 
resolution of the high altitude structure of the disk, if serious 
errors in the cooling rates derived from radiative emission 
from the disk photosphere are also to be avoided. 

At the orbit radii most relevant for planet formation 
(<20 AU), disks will be optically thick in the sense that 
the optical depth, r, calculated from infinite distance to the 
disk midplane is large. This is important because it means 
that the cooling rate of packet of disk material will be well 
modeled by a blackbody cooling law whose temperature is 
defined at the disk photosphere. If the altitude of the photo- 
sphere is known only approximately, due to low resolution or 
simple miscalculation, the cooling rate will be similarly ap- 
proximate because the temperature structure near the pho- 
tosphere is not known precisely. Comparatively small errors 
in temperature will lead to much larger errors in the local 
cooling rate because of the T 4 proportionality of the black- 
body cooling function. 

Accurate knowledge of the cooling rate and its asso- 
ciated photosphere temperature is important because the 
dynamical balance of heating and cooling processes in disks 
is quite precarious: most of the heating and cooling sources 
are capable of removing or replacing essentially all of the 
disk' s thermal energy ov er the course of only a few or- 
bits iDurisen et al. \\200db- S pecifi cally for the question of 
disk fragmentation. iGammiel ll200 ill shows that cooling rates 
faster than ~ 3/fi, where f2 is the local orbit frequency, lead 
to disk fragmentation, but longer cooling rates may not in 
a loca l calculation. Later global calculations of R ice et al. I 
(2003) confirmed this result for low mass disks, where they 
found that a cooling time of 5/57 was long enough to prevent 
fragmentation, but a cooling time of 3/fi was not. Simula- 
tions of higher mass disks with cooling rates of 10 and 5 
fi _1 showed similar changes in behavior. From these results 
it is clear that an error in the cooling rate as large as a 
factor of two will be sufficient to suppress or enhance frag- 
mentation in a simulation which would otherwise follow a 
very different evolutionary path. A error of factor of two in 
the cooling rate will be equivalent to an error of ~ 20% in 
the photosphere temperature, due to the dependence 
of the cooling on the temperature. Since an error of such 
magnitude may seriously alter the balance between cooling 
and heating, a more restrictive limit must be set to ensure 
that any fragmentation that does occur is not subject to nu- 
merical error. If we set an arbitrary standard that a cooling 
rate accurate to 25% will be sufficient to ensure results are 
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Figure 24. Vertical density (top), temperature (middle) and 
optical depth (bottom) profiles of a disk at 10 AU, with local sur- 
face density of 500 gm/cm 2 and a midpl ane temperature of 50 K, 
using the model of iNelson et al. 1 120001) to determine structure. 
The altitude coordinate is normalized to the value of the isother- 
mal scale height H lso determined at the midplane. Dotted and 
long dashed vertical lines on each curve are placed to straddle 
the photosphere surface at r = 2/3 (horizontal dashed line), at 
distances of -Hj so /10 and Hj so /20, respectively. Over these dis- 
tance, the temperature deviates ~ 15% and ~ 7% above or below 
the true photosphere temperature. 

not contaminated by numerical errors, then we will in turn 
require a photosphere temperature accurate to ~ 5 — 6%. 

Given the requirement as stated, it remains to convert 
the constraint on the photosphere temperature determina- 
tion to a constraint on required spatial resolution of the ver- 
tical disk structure. One approximate model of the vertical 
struct ure of a disk is available from the work of lNelson et al. I 
(2000), in which an isentropic vertical structure is assumed, 
for a locally plane parallel, self gravitating disk. Figure 1241 
shows temperature and optical depth profiles as a function 
of altitude for a region of a disk expected t o be interest- 
ing fo r disk fragmentation, derived from the lNelson et al~\ 
2000 model. Several points are of interest in evaluating the 
figure. First, although the density structure follows an ap- 
proximately Gaussian structure, at least in coarse outline, 
it is dramatically compressed relative to that expected for 
an isothermal disk. No material is present above an altitude 
of ~ 1.7-ffiso, due to the different physical assumptions em- 
ployed in deriving the mass distribution. Second, this con- 
figuration is optically thick, with an optical depth to the 
midplane of ~ 100, meaning that a disk photosphere surface 
may in fact be defined. 
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For the purposes of evaluating the radiative energy 
losses, a third point will be of most importance. Namely, 
that near the disk photosphere surface both the density and 
temperature are changing rapidly with altitude. A calcu- 
lation that resolves these profiles too coarsely will therefore 
also only coarsely resolve the photosphere surface itself, pos- 
sibly exposing the hot disk interior to space. For example, if 
the particle density near the photosphere for an SPH simu- 
lation is such that smoothing lengths are h m Hi so /10, then 
to a first approximation, the photosphere surface itself will 
be resolved on the same spatial scale because interparticle 
separations are themselves ~ h. A particle actually located 
anywhere within h of the actual photosphere surface may 
be determined to define the photosphere surface, simply be- 
cause no other particle happened to be at a higher altitude 
there. An identical argument holds for grid based simula- 
tions, with the size of a zone at the photosphere surface 
replacing smoothing length. 

Fieure l^H shows that uncertainties in the altitude of the 
photosphere surface of ~ Hi so /10 results in an increase or 
decrease of the derived photosphere temperature of ~ 15%. 
Doubling the linear resolution to ~ _ffi so /20 decreases the 
uncertainties to ~ 7%, comparable to our requirement of 
errors no larger than ~ —5 — 6%. We therefore conclude 
that simulations including radiative cooling must resolve the 
disk photosphere surface at a spatial scale no coarser than 
~ Hiso/20, when structure models similar to ours are em- 
ployed. 

After completing the exercise and deriving this require- 
ment, we must also immediately point out that the quoted 
required resolution is not universally applicable without 
modification. Any constraint on spatial resolution in the con- 
text of radiation will be considerably complicated by the fact 
that the vertical temperature profile will be significantly in- 
fluenced bothjjvjjrrferential dissipation of shocks at high al- 
titudes l|Pickett et al. 2003), by radiative heating by exter - 
nal sources like the central star JChiang fc Goldreichlll997l) 
and, if any remains at the ti me of the simulation ; the sur- 
rounding molecular cloud. As lChiang fc Ooldreichll997l also 
note, wavelength dependent opacities will also play a role, 
to the extent that long wavelength radiation is able to cool 
the disk interior effectively, even while shorter wavelength 
radiation remains trapped. 

In spite of these cautions, we remain convinced that the 
above exercise remains a valuable illustration. It is based 
only on the condition that the transition between optically 
thick and optically thin material is well resolved and, al- 
though it is applied to a particular model of the vertical 
structure, it is not tied to it specifically. Given the overall 
gross vertical structure of disks, in which densities are high 
at the midplane and drop steeply as altitudes increase, op- 
tical depth profiles similar to that in figure |2"H will occur in 
other models as well. Because the optical depth is a very 
steeply falling function of increasing altitude, the transition 
region will be narrow, as in the model above. Whether that 
narrowness can be translated directly into a constraint on 
the spatial resolution will depend on the details of the tem- 
perature profile there, and its influence on the accuracy of 
radiative energy transfer. The specific criterion as applied 
above presumes a temperature profile decreasing with alti- 
tude as the density decreases, T oc p 7_1 . If instead a much 



flatter profile exists across the transition, a less restrictive 
criterion may still permit sufficiently accurate evolution. 



6 CONCLUDING REMARKS 

After exploring the conditions required for simulations to 
produce numerically valid results, we turn now to summa- 
rizing those conditions, to comparing the work done here 
with previous work in the literature, and to commenting on 
additional requirements beyond numerics that are required 
for simulations to accurately reproduce the evolution of real 
systems. 

6.1 Summary of our main results 

Our first result is to extend the numerical criterion for 
the validity of simulations involving collapse in clouds and 
based on the analyt ic Jeans collapse formalism (BB97, 
iTruelove et al. I lll997f) ). to the case of disk systems where 
the geometry of the mass distribution is both flattened 
and rotating so that the Jeans treatment is inapplicable. 
Whereas a minimum spatial resolution per Jeans wavelength 
or minimum number of particles per Jeans mass is required 
in the cloud collapse case for grid or particle based simula- 
tions, respectively, in simulations of disks we require a min- 
imum spatial resolution per 'Toomre wavelength' or mini- 
mum number of particles per 'Toomre mass', for grid based 
or particle base simulations, respectively. The two forms of 
the new criterion are given in equations l5l and 1111 

We also note that simulations failing to resolve one or 
the other condition yield quite different outcomes. In the 
case of the Jeans condition, BB97 conclude that insufficient 
resolution will delay or suppress numerically induced frag- 
mentation with softening and smoothing equal, and only en- 
hance it when smoothing is much larger than softening. Our 
results demonstrate that a simulation that fails to resolve the 
Toomre condition is characterized by enhanced fragmenta- 
tion compared to one that is well resolved, even when soft- 
ening and smoothing are equal. This fact allows an empirical 
determination of the minimum resolution required in a sim- 
ulation in order to ensure fragmentation that develops is not 
of numerical origin. For 2D SPH simulations, we have deter- 
mined the minimum number of particles per Toomre mass 
as six times the average number of neighbors (~ 20) used for 
the realization of the hydrodynamic quantities themselves. 

We have not determined a value of N ICSO for 3D mod- 
els, but its seems unlikely to be dramatically different from 
2D, since the same physical system is being modeled with 
the same numerical technique. To the extent that differences 
may be present, we expect that the criterion will be more 
conservative (more particles required) based on two argu- 
ments. First, unstable waves in the Toomre analysis are in- 
trinsically two dimensional, so resolving their third dimen- 
sion is meaningless: small particle separations or small grid 
spacing in that coordinate do not contribute to increasing 
the resolution of the waves themselves. Using the same em- 
pirical fragmentation/non- fragmentation limiting condition 
used in section l3~3*l to establish a required particle count, will 
therefore overestimate the required resolution because some 
neighbors will be counted only for purposes of establishing 
a value for N le so, but will not actually serve to increase 
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the resolution insofar as the Toomre condition is concerned. 
Secondly, although the Toomre criterion implicitly assumes 
validity of the hydrodynamic equations, in practice the nu- 
merical method may fail to do so (sectional , with the conse- 
quence that pressure forces are underestimated. Although a 
vertical resolution requirement is nominally an entirely dif- 
ferent condition from resolving instability wavelengths, we 
see no practical manner in which they might be separated 
from each other in the empirical tests we describe. Therefore, 
calibration of the requirements for the Toomre condition in 
3D will implicitly incorporate both requirements from the 
wave analysis and from the vertical resolution condition, 
yielding a resolution requirement defined by contributions 
from both sources. 

A second important result of our work is to reconcile a 
variety of mutually inconsistent definitions of the Jean col- 
lapse stability criteria made by different authors, due to the 
use of different definitions of the Jeans mass itself. We show 
that using the same definition of the Jeans mass, a particle 
based simulation will require approximately 600 particles to 
resolve one Jeans mass, while a grid based simulation will re- 
quire about 64 zones (i.e. 1/ J 3 where J is defined in equation 
and J ~ 1/4). The minimum resolution of 600 particles 
per Jeans mass corresponds to approximately 12 times the 
average number of neighbor particles (~ 50) used by BB97 
in their SPH code. This value is a factor of two higher than 
the factor of six times the average number of neighbors as 
derived from the 2D disk simulations explored in section l3~2l 

Third, we discuss the importance of the implementa- 
tion of interparticle forces in particle based simulations of 
self gravitating disks. We show that a naive translation of 
the 3D method for gravitational softening in SPH into 2D 
leads to an unphysical, finite gravitational force at zero sep- 
arations between particles. In consequence, significant par- 
ticle 'pairing' occurs as simulations proceed, resulting in an 
effectively time dependent resolution. We consider two alter- 
natives to remedy the behavior. First, we artificially mod- 
ify the derivative of the SPH kernel according to the TC92 
prescription, generating a finite, repulsive pressure force at 
zero separation which compensates for the attractive grav- 
itational force. Second, we artificially modify the gravita- 
tional softening itself, to produce instead a zero force at 
zero separation, as is physically appropriate. 

Both alternatives are effective in reducing the likelihood 
of disk simulations to produce fragmentation and of parti- 
cles in those simulations to become paired. Although the 
TC92 variant appears to be slightly more effective (see e.g. 
section 14.311 . we view it as marginally the less preferable of 
the alternatives because it too produces a finite force at zero 
separation, albeit of opposite sign so that a sum of gravita- 
tional and pressure forces yields a small net force. Finite 
forces at zero separation, whether or not they are combined 
with others to produce a near zero sum, must be considered 
a less physical and therefore less desirable alternative to a 
solution which retains the collisionless nature of the model 
more directly. 

Fourth, we show that force imbalances on the spa- 
tial scale of the smoothing and softening develop when the 
two length scales are not equal. While acting only at very 
small spatial scales, these imbalances can induce large scale 
changes in the outcome of simulations where they exist, in- 
cluding artificially induced fragmentation in otherwise stable 



systems. To avoid such induced outcomes, we conclude that 
the gravitational softening length and the hydrodynamic 
smoothing length of each particle must be very nearly iden- 
tical. As noted above, our conclusion for disk simulations is 
similar to that made by BB97 in their study of collapsing 
clouds, but stronger in the sense that it holds even when the 
systems are, and are expected to remain, properly resolved 
according to the Toomre criterion. 

In neither the disk or cloud case however, does our con- 
clusion come without price. Specifically, in order to better 
resolve the flow in both high and lower density regions of the 
same simulation, most modern SPH codes set the smoothing 
length dynamically as a function of the local flow. Setting 
the smoothing and softening lengths equal means that the 
softening length also must vary in time. This is important 
because allowing temporally varying softening is equivalent 
to explicitly allowing a violation of conservation of energy in 
the simulation. While we believe the magnitude of violation 
to be small, its effect must be considered when evaluating 
self gravitating particle based simulations. 

Finally, we discuss the importance of adequately resolv- 
ing the vertical structure of 3D disk simulations. We find 
that insufficient resolution of the vertical structure leads to 
substantial deficits in the realized midplane densities com- 
pared to the predictions of analytical models. We find that 
resolution comparable to at least ~ 4 smoothing lengths 
per disk scale height in the disk midplane are required to 
accurately fit the coefficients in analytic formulae for the 
vertical structure of non self gravitating isothermal disks in 
particle simulations. Restated in terms of the full thickness 
of the disk, the criterion is equivalent to the statement that 
at least ~ 17—20 particles are required per vertical column. 
While we have not studied the effects of resolution in 3D 
grid based simulations, we expect similar criterion to hold 
for them as well. 

The consequences of failures to resolve the vertical 
structure in self gravitating systems will be that pressure 
forces will be incorrectly determined. Any balance between 
them and self gravitation will therefore be falsely biased in 
favor of gravitational fragmentation. Presumably, the bias 
will distort the evolution of waves in the disk (whether self 
gravitating or not) even when fragmentation does not occur, 
because the density and pressure deviations imply a loss of 
fidelity of the solution to the equations of hydrodynamics 
provided by the numerical method. 

We discuss the fact that while satisfying the criterion 
above is necessary for simulations to accurately evolve cir- 
cumstellar disks, it may not be sufficient for them to do so 
alone, depending on the goals of the researcher. Simulations 
for which the criterion is satisfied everywhere may still under 
resolve high altitude structure of the disk, a fact which will 
be of great importance for models in which radiative trans- 
port is included. Such models must also resolve the temper- 
ature structure at the disk photosphere in order not to be 
contaminated by unacceptably large numerical errors. We 
provide an example of the requirements when an isentropic 
vertical structure is present, for which we derive the com- 
putationally extremely demanding condition that the pho- 
tosphere must be resolved at a spatial scale smaller than 
~ H/20 to avoid serious errors. We point out however that 
the exact requirement will be highly sensitive to the details 
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of the model and its effects on the temperature, density and 
opacity structures at high altitudes. 

6.2 Comparisons to particle based simulations in 
the literature 

Our test prob lem clearly shows tha t the fragmentation pro- 
cesses seen m iNelson et all dl99Sl) were due to insufficient 
resolution, and we believe that many other works in the lit- 
erature both before and since may be similarly affected. Of 
the work dis cussing particle based sim ulations, the simula- 
tions done bv lMaver et al. I ([2002 , 2004) are among the high- 
est resolution simulations so far published, evolving either 
2 x fO 5 or fO 6 particles (in various runs) in three dimen- 
sions. It may therefore be particularly useful to contrast our 
test problems with the conditions employed by them, their 
results and conclusions. 

They perform a series of SPH simulations with several 
disk masses between M D /M„=0.075 and M D /M*=0.125, 
each with minimum Toomre Q values between 0.8 and 1.9. 
They show that the more massive disks fragment on a time 
scale of a few hundred years (about 10 orbits in the region 
where the clumps form) when the minimum Q values in the 
disk are below Q m i n ~ 1-4. They also find that clump forma- 
tion is more vigorous with higher resolution simulations of 
the same initial condition that with lower. An isothermal gas 
equation of state was used initially for most of the models, 
and after a critical over density was exceeded, was switched 
over a small transition time to an adiabatic evolution that 
included heating due to compression, viscosity and shocks. 
Other models used an adiabatic equation of state for all 
times, and resulted in fragmentation that was much weaker 
or non-existent. 

We first attempt to obtain some indications of the de- 
gree to which the flow is well resolved by the condition in 
equation 1111 at different times during the evolution, using 
their discussion of the physical features of the flow and res- 
olution. For example, prior to clump formation the most 
unstable wavelength is quoted as ~ 2 AU, including a cor- 
rection factor for 3D effects. The flow is clearly well resolved 
by the simulation at this early stage, but the value of the 
Toomre wavelength 4 at this time is only representative of 
the mildly non-linear regime of the background flow. It is 
not indicative of the later, more strongly non-linear regime 
near the forming clump, for which both the locally applied 
Jeans based or Toomre based resolution criteria are most 
useful. 

Well after a clump has formed, when the clump has 
become a well defined, condensed and independent en- 
tity rather than a part of the disk, we expect that the 
Toomre based crite rion will become invalid. At this time, 
iMaver et al\ [2004 state that the resolution employed in 
their simulations was sufficient to resolve the Jeans mass (for 
which the underlying analysis is valid if the clump is much 
smaller than a disk scale height) with a few hundred to a few 
thousand particles, even though the densities there reached 
as much as ~ 10 6 higher than the background. An important 

4 Note that their definition of the 'Toomre wavelength' is that 
of equationl^l but that their later definition of the most unstable 
wavelength matches that of equation 151 



detail of the implementation o f lMav er et al. 120041 is that the 
number of neighbors used for realizing the hydrodynamic 
quantities is fixed at 32, a value common throughout the 
cosmological SPH literature, but significantly smaller than 
the value of ~ 50 used (but not specifically recommended) 
by BB97. Correcting for this systematic difference would not 
appear to alter the statement that the Jeans mass is prop- 
erly resolved however, because while the number of particles 
per Jeans mass would be lower, it would still remain above 
the ~ 100 quoted by BB97. 

We can make relatively few direct inferences about in- 
termediate times when the clump has begun to form but 
is still unbound or only weakly bound. At such times the 
local value of Q may fall below one half, and we expect a 
much stronger resolution constraint from the Toomre crite- 
rion than from the Jeans criterion (see equation 1161 . espe- 
cially if the approximation that fl « re is not made. There- 
fore, while the clump may be well resolved by the Jeans cri- 
terion, it may not be similarly well resolved by the Toomre 
criterion. Based solely on the number of particles and the 
similarity of the initial conditions to our own test problem, 
we conclud e that the high reso lution (1 million particle) sim- 
ulations of ilVIaver et al. 1 l2004l appear to be well resolved by 
our Toomre based resolution criterion. We can make a simi- 
lar, but much more tentative conclusion about the lower res- 
olution (2 x 10 5 particle) versions in the same work. A more 
detailed examination of the simulations at both resolutions 
in the context of the Toomre resolution criterion would be 
desirable, particularly in the context of the sensitive depen- 
dence of the resolution requirement to the initial minimum 
value of Q and on the local value of Q in the fragmenting 
regi on. 

IMaver et al. \\20o4 advocate a fixed gravitational soft- 
ening whose length scale is half of the initial hydrodynamic 
smoothing length of the particles. With this ratio of soft- 
ening to smoothing, BB97 show that the gravitational force 
can exceed the pressure force by a factor of as much as seven 
on size scales of h/2, and by a factor of 2-3 on scales of h. 
The actual length scale used is determined (Mayer 2004 pri- 
vate communication) as an average of the smoothing lengths 
in the 'flat' temperature region of the disk, corresponding 
approximately to the outer half of its radial extent, mean- 
ing that the outer half of their disks will have small scale 
force imbalances that originate in the numerical treatment 
of softening and smoothing. 

The results of our simulations show the consequences 
of such imbalances in disk simulations. In regions where the 
hydrodynamic smoothing length exceeded the gravitational 
softening length, fragmentation was strongly enhanced over 
that of the variable softening case, even when sufficient nu- 
merical resolution was used according to the Toomre crite- 
rion. Since the force imbalance is not due to any physical 
origin, but only to a breakdown in the simulation method 
itself, we conclude that clumping observed in the simulations 
may not be an accurate depiction of an underlying physical 
system, and should be reexamined to ensure that they are 
not subject to this numerical flaw. 

Of still greater concern is the resolution of their simu- 
lations in the context of the discussion i n section 15.61 Both 
the morphology of the system modeled bv lMaver et ~al. 120041 
and the physics itself are different than those employed by 
our test models, however it seems unlikely that these modifi- 
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cations will alter the required resolution by the large factors 
needed to conclude that their vertical structure is well re- 
solved, especially for their most unstable (lowest Q) models. 
The extent to which any violation of the criterion actually 
affects fragmentation is unclear however, ft is likely that a 
definitive statement of the effects will require a complete 
study of the problem at extremely high resolution, making 
the question quite costly to address at present. 

Several other recent papers have discussed fragmenta- 
tion in disks. Similar concerns regarding the vertical reso- 
lution will apply to many of them. For example, in iBosj 
(2004), as well as his earlier, similar work, the vertical co- 
ordinate (the spherical coordinate 6 in his case) is resolved 
with a total of 23 zones. Although they are asymmetrically 
arranged to concentrate resolution where they are needed 
most, only some ~10-12 lie within ten degrees {H/R m 0.17) 
of the midplane, where most disk matter will be located. Al- 
though our vertical resolution requirement has been tested 
only with SPH simulations, there is no reason to believe that 
an analogous requirement will be required for grid based 
simulations as well, since they must solve the same set of 
equations. Naively applying the ~ 4 smoothing lengths per 
scale height as a ~ 4 grid cells per scale height criterion di- 
rectly to the Boss simulations, we observe that the vertical 
structure may be very near the required limits. The exact 
number of cells will undoubtedly be different than a straight- 
forward carryover however, so the extent to which Boss's or 
any other any grid based simulations will be affected by the 
erro rs arising from def icient vertical resolution is unclear. 

iRice et al. I i2003f) showed that for sufficiently fast cool- 
ing rates, clumping could occur in disks resolved with 
2.5 x 10 5 particles. Due to their cooling prescription, disks in 
their work dropped to globally averaged Q values near unity, 
well below those in our simulations, for which we determined 
erroneously low midplane densities developed. In later work 
discussing similar systems (with overall morphologies iden- 
tical to those i n our study, but with additional physics), 
lLodato &z Ricei i2004h claim that their vertical structure is 
resolved with ~ 5 smoothing lengths per scale height, ft is 
difficult to reconcile this result with our own figure l2~2l es- 
pecially given the much lower Q values (see th eir figure 3) 
realiz ed in their study and in the earlier work of IRice et al. I 
120031 . It appears likely however, that low vertical resolution 
has likely influenced their results through the artificially low 
pressure forces induced by density deficits. 

Apart from the issue of vertical resolution, we note that 
their overall resolution (in 3D) is similar to that in our high- 
est resolution 2D test problem. The authors have not dis- 
cussed resolution issues in the context of either Jeans or 
Toomre stability, so we can draw few inferences about the 
viability of their results under either the criteria of equation 
Hlor llll In the context of the latter condition however, the Q 
values throughout their simulations rapidly approach values 
near unity. This is important because equation ll2l shows that 
our resolution requirement is quite sensitive to the value of 
Q, so that much higher resolution will be required for later 
evolution even if the initial condition is known to be well 
resolved by the Toomre co ndition. 

Other simulations of Lufkin et al. I (l2004h perform a 
large parameter study of disks already containing one planet, 
in order to explore the possibility that its present may trig- 
ger the formation of further objects. For their parameter 



study, they use simulations with 10 5 particles, but with vary- 
ing minimum Q or initial planet mass. The initial conditions 
of these disks were in most respects quite similar to those 
of Mayer et al. , however the resolution employed was far 
lower in order to facilitate the large parameter study. Given 
the much lower resolution, it appears likely that these au- 
thors have under resolved the vertical structure of the disk, 
resulting in incorrectly low densities and pressures and the 
possibility of artificially enhanced fragmentation. While an 
attempt to duplicate their study is beyond the scope of this 
work, we have run selected simulations sampling their pa- 
rameter space and have observed that they resolve (or fail 
to resolve) the Toomre condition at similar points in the pa- 
rameter space for which they conclude that clumping was 
induced or not induced by the passage of the seed planet. 
Due to the probable violation of both criteria, we therefore 
believe that their conclusions have been contaminated by 
insufficient resolution and must be reexamined. 

6.3 Beyond numerics: what is required for the 
viability of the disk instability model for 
Jovian planet formation to be either verified 
or falsified 

In this work, it has not been our purpose to discuss all 
of the requirements for simulations modeling gravitational 
fragmentation in disks to produce correct and physically rel- 
evant results. Such a description of a real physical system 
requires not only numerically valid simulations, but also rel- 
evant initial conditions and a correct and complete physical 
model. Numerically valid simulations will be interesting only 
to the extent that they model real systems with real physics. 

Of particular interest for the physical understanding 
of disk evolution is the fact that our highest resolution 
SPH test model did not produce any clumps although it 
was evolved for much longer in time than its lower resolu- 
tion cousins. This is an important physical result because 
the test problem employed a locally isothermal equation of 
st ate, which is usually (and erroneously-see the discussion 
in lNelson et al. I i200flh : iPickett et al. I ll2003l) ) equated with 
the strong cooling limit where fragmentation is nearly in- 
evitable. The fact that clump formation did not occur in 
this limit, in a properly resolved simulation, may place the 
disk instability model for planet and brown dwarf formation 
on considerably less solid ground. 

At the same time, the Boss simulation, implementing 
conditions that are much less biased toward clumping, did 
produce clumps while remaining numerically valid accord- 
ing to the Toomre condition, though some doubt may re- 
main regarding the resolution of the vertical structure. His 
simulations include a description of radiative cooling pro- 
cesses in 3D, and he concludes that the clumping is due 
both to the efficient radiative cooling and to the fact that 
convection in the disk's vertical direction is efficient. Con- 
vection is able to transport thermal energy out of the opti- 
cally thick disk midplane to its photosphere surface at higher 
altitudes, so that even optically thick regions like forming 
cl umps cool effic i ently. On the other hand, previous 2D work 
o flNelson et al. I feOOOTl , which assumed efficient vertical con- 
vection as a consequence of the radiative cooling prescrip- 
tion employed, did not produce clumps. Moreover, the to- 
tal radiated energy emitted by the disks modeled in those 
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simulations seriously underestimated the values observed in 
real systems because they were too cold. In a more com- 
plete physical model of the disk energy budget of heating 
and cooling processes, the disks would have to be warmer in 
order to reproduce the observations, making clumping even 
less likely. 

Other fully 3D radiative transfer (iMeiia et al. 11200.1 ) 
models confirm neither the Boss conclusion about the effi- 
ciency of vertical convection nor the clumping that results 
from it. The authors continue to investigate the origin of 
the discrepancies, but no certainty yet exists. Possibilities 
include both differences in the numerical treatments, for ex- 
ample of the boundaries, the resolution of both the Poisson 
solver and the hydrodynamic codes themselves, and of the 
physical models, for example the treatment of external ra- 
diation field, what opacities are used in the radiative trans- 
port and the equation of state for the gas. The important 
point to take from these discrepant results is that not only 
do questions of numerical validity remain to be addressed, 
but also the physical models. Although progress has cer- 
tainly been made, to date an insufficient understanding of 
what physical processes are important for fragmentation has 
been developed. Determination of the true viability of the 
disk instability model for Jovian planet formation therefore 
still lies in the future. 
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